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1.0  SUMMARY 


Road  maps  can  provide  information  on  not  only  topographical  conditions  but  also  on  how 
materials  and  people  move.  Historically,  available  information  pertaining  to  road  maps  consisted 
of  static  images  of  preordained  routes.  Even  now,  Google  maps  and  hand  held  Global  Positioning 
Systems  (GPS)  represent  a  static  view  of  road  networks,  requiring  either  recapturing  images  or 
manually  updating  units.  However,  in  order  to  have  more  current,  information  rich  representations 
of  transportation  networks  or  road  maps,  the  use  of  kinematic  information  provided  by  sensors 
such  as  Ground  Moving  Target  Indicator  (GMTI)  data  is  required.  The  data  these  GMTI  sensors 
supply  is  not  only  a  static  representation  of  a  target  but  also  the  kinematics.  The  approach 
employed  for  synthesizing  the  received  data  into  a  complete  estimate  of  the  road  network  is 
through  the  use  of  the  Hough  Transform,  to  identify  line  segments  which  collectively  represent 
the  road  network.  The  Total  Least  Squares  is  used  to  characterize  the  uncertainty  associated  with 
this  representation  as  well  as  provide  a  more  accurate  estimate.  The  uncertainty  in  each  of  the 
estimated  parameters  (i.e.  slope  and  intercept  of  a  line)  can  then  be  approximated  by  the  Cramer 
Rao  lower  bounds.  Finally  the  identified  segments  are  merged  and  connected  to  provide  a  more 
complete  representation  of  the  road  network.  The  approach  used  here  is  iterative,  for  example, 
when  new  data  within  the  area  of  interest  is  received,  a  better  estimate  of  the  road  segment  can 
be  obtained  and  the  overall  road  network  is  updated  and  allowed  to  grow  in  all  directions. 
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2,0  INTRODUCTION 


The  need  to  develop  aeeurate  estimates  for  road  networks  is  important  in  both  military 
eivilian  applieations.  Synthetic  Aperture  Radar  (SAR)  and  Ground  Moving  Target  Indieator 
(GMTI)  data  is  often  proeessed  and  analyzed  to  produee  sueh  networks.  SAR  produees  images 
of  varying  intensity  whieh  ean  be  proeessed  to  separate  buildings,  roads,  and  terrain.  However, 
SAR  is  only  able  to  deteet  prominent  existing  features  [1].  In  other  words,  SAR  will  only  deteet 
a  road  if  there  is  a  distinet  outline  of  sueh  a  path.  GMTI  on  the  other  hand  traeks  moving  targets 
and  relays  the  latitude  and  longitude  eoordinates  as  well  as  the  veloeity  in  both  direetions.  The 
disadvantage  of  GMTI  however,  is  the  neeessity  of  a  moving  target,  should  the  target  stop  or 
be  obstrueted  in  any  way  from  the  sensors,  the  traeker  will  lose  the  target  for  the  duration  of 
the  obstruetion  [2].  Many  of  the  eurrently  available  algorithms  rely  on  information  from  pre¬ 
existing  road  maps[  [1],  [3],  [5],  [7]],  however  in  many  seenarios  the  availability  of  this  a-priori 
information  is  limited  and  inaeeurate.  In  additional  situations  there  is  no  existenee  of  road  maps, 
sueh  as  in  times  of  eonfliet  in  desert  regions.  Therefore  the  need  for  an  algorithm  whieh  ean 
aeeurately  estimate  road  networks  in  a  timely  manner  is  of  great  importanee.  Furthermore,  there 
is  a  laek  of  a  quantifiable  measure  of  the  aeeuraey  of  the  extraeted  road  estimates.  Several 
available  algorithms  use  a  ’Completeness”  and  ’’eorreetness”  measure  whieh  is  a  eomparison  of 
the  extraeted  road  network  and  the  aetual  network  [  [7],  [8],  [9],  [11]],  however  as  previously 
stated  in  many  situations  there  is  no  available  truth  so  these  methods  are  not  relevant. 

2.1  Related  Work 

Koeh,  Roller  and  Ulmke  utilize  Multiple-Hypothesis  Traeking  (MHT)  to  develop  the  estimates 
of  the  target  state  veetor  Xk  whieh  is  eomposed  of  the  ground  eoordinates,  rk,  and  their 
eomponent  veloeities,  f^,  at  time  k  [2].  In  most  seenarios  the  z-eomponent  of  both  the  ground 
eoordinates  and  the  veloeity  are  ignored  or  set  to  zero.  The  MHT  algorithm  is  outlined  in  Figure 
1. 

MHT  is  a  stoehastie  logie  deeision  maker  based  on  both  eurrent  and  future  observations.  Suppose 
there  exist  two  traeks,  and  at  the  time  step,  k,  three  observations  are  reeeived.  Multiple  hypotheses 
are  made  as  to  whieh  traeks  these  observations  eorrespond  to,  either  one  traek,  both  traeks,  or 
neither  of  the  traeks  are  plausible  hypotheses.  Eaeh  of  these  hypothesis  are  then  propagated  into 
the  future,  and  onee  additional  observations  are  obtained  these  hypotheses  are  evaluated  and  the 
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USER 


Fig.  1:  MHT  Algorithm  [3]  pp.  7 


most  probable  track  is  accepted.  The  first  box  in  Figure  1,  termed  ’’Gating”,  forms  a  region 
based  on  the  covariance  matrices  of  the  predictions  and  only  hypotheses  within  this  region  are 
accepted,  this  reduces  the  initial  number  of  hypotheses.  For  the  MHT  algorithm  the  target  state 
is  modeled  by  a  linear  Markov  process  which  is  defined  such  that  the  state  at  time  tk+i  is 
dependent  on  only  the  preceding  time  step,  tk-  The  measurement  function  is  also  assumed  to  be 
linear  in  terms  of  the  state.  Retro-diction  allows  for  the  smoothing  of  the  hypotheses,  it  does 
this  by  recalculating  the  previous  estimates  with  all  available  measurements.  In  this  algorithm  all 
roads  are  modeled  by  linear  segments,  additional  segments  are  added  depending  on  the  linearity 
of  the  road.  In  order  to  describe  the  linearity  of  the  road,  the  directional  vector  of  the  velocity, 
ri\k  and  the  difference  in  the  ground  coordinates,  ri+i|k  -  ri|k  are  considered.  Where  I  denotes 
the  node  and  k  is  the  time  step.  Figure  2  depicts  a  graphical  representation  of  a  curved  road 
segment  where  additional  nodes  must  be  incorporated. 

The  angle  of  the  distance  vector,  ip^k  and  the  angle  made  by  the  tangential  velocity  vector  and 
the  distance  vector,  denoted  by,  (f)i\k  are  used  to  develop  Equation  1,  the  decision  criteria  for 
including  additional  nodes  for  line  segments. 


$ 


>  K 


l\k 


(1) 


Where  k  is  the  decision  parameter  and  after  experimentation  is  selected  to  be  one  and  ^i\k 
denotes  the  covariance  matrix  of  the  velocity. 

Baumgartner,  Hinz,and  Wiedemann  developed  a  semi-automatic  method  for  line  extraction 
from  images  [4].  The  user  is  required  to  identify  an  initial  segment  which  the  tracker  will 
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Fig.  2:  Linearity  Determination  [2]  pp.  203 


use  to  then  grow  the  entire  network  and  should  the  automatic  tracker  encounter  an  error,  the 
user  is  required  to  take  action.  Two  main  classes  of  problems  can  be  encountered  during  the 
implementation  of  the  road  tracking  algorithm.  The  first  case  is  when  the  tracker  reaches  the 
edge  of  the  image.  The  second  case  is  when  the  internal  confidence  of  the  tracker  is  below  a 
desired  level.  This  algorithm  works  well  for  distinct  road  like  features  within  the  image.  There 
are  five  steps  performed  during  the  course  of  the  estimation  of  the  road  network  defined  below: 
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1)  Generate  a  reference  profile 

2)  Predict  the  next  position  of  the  road  axis 

3)  Estimate  the  position  of  the  next  point  on  the  road 

4)  Check  for  stopping  criteria,  if  none  return  to  (2)  otherwise  (5) 

5)  Request  user  interaction 

This  process  is  repeated  until  all  desired  roads  have  been  identified. 

Zhou,  Venkateswar,and  Chellappa  begin  the  extraction  of  roads  with  an  edge  detector  to 
identify  the  boundaries  of  such  roads  [5].  This  produces  edge  pixels  which  are  then  utilized  to 
identify  linear  features  in  an  image.  The  linear  feature  extraction  is  broken  up  into  three  steps. 
First  the  gradient  direction  of  the  edge  pixels  must  be  determined.  For  simplification  purposes,  the 
gradient  direction  is  discretized  into  eight  evenly  spaced  values  between  -180  and  180  degrees. 
To  determine  whether  two  edge  pixels  are  collinear  the  template  of  the  edge  pixel  is  compared 
against  a  database  of  available  templates.  The  templates  for  four  of  the  eight  gradient  directions 
are  shown  in  Figure  3. 


direction; 1 
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Fig.  3:  Pixel  Templates  [5]  pp.  91 


Each  edge  pixel  is  then  separated  into  the  corresponding  line  structure.  The  second  step  is  to 
connect  the  identified  lines.  A  small  neighborhood  is  examined  around  each  of  the  endpoints  of 
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a  line.  Two  criteria  must  be  met  in  order  for  the  lines  to  be  merged.  First  the  direction  of  the 
second  line  must  be  similar  to  the  first,  and  second,  the  magnitude  of  the  second  line  must  lie 
between  half  and  twice  the  magnitude  of  the  first  line.  This  step  is  performed  twice,  once  for  a 
direction  delta  of  8  degrees  and  second  for  a  delta  of  13  degrees.  The  third  and  final  step  in  the 
extraction  is  the  extension  of  line  segments  to  form  comers.  This  is  easily  done  by  determining 
if  two  lines  will  intersect  due  to  their  directional  gradients  and  if  so,  they  are  extended  until 
they  intersect.  This  corner  extraction  method  however  is  only  viable  for  two  intersecting  lines, 
should  there  be  more,  the  comer  will  fail  to  be  an  accurate  estimate. 

Amo,  Martinez,  and  Torre  extract  roads  from  aerial  images  using  a  region  competition  al¬ 
gorithm  [6].  This  method  is  semi-automatic  in  which  the  user  selects  points  along  the  road. 
The  amount  of  points  required  varies  directly  with  the  curvature  of  the  road  segment.  These 
user  selected  points  must  be  located  at  large  curvature  changes  and  also  where  the  intensity  of 
the  image  classifies  a  region  of  the  road.  The  user  must  also  initialize  the  maximum  allowable 
distance  between  two  points,  S,  and  the  minimum  allowable  angle,  0,  between  the  triplets  of 
points.  These  user  selected  points  are  then  interpolated,  and  where  needed  additional  points  are 
automatically  added,  and  a  B-spline  is  produced,  corresponding  to  the  centerline  of  the  road. 
With  this  centerline  the  outer  edges  of  the  road  are  determined  by  simply  copying  the  spline 
and  moving  it  until  the  statistical  properties  of  the  region  no  longer  coincide  with  the  statistical 
properties  of  the  road,  designating  a  relative  bound  of  the  road.  These  newly  produced  edge- 
splines  are  however  inaccurate  and  the  use  of  region  competition  refines  these  initial  estimates 
since  widths  vary  and  road  edge  paths  may  vary  as  well.  These  so  called  regions  are  homogeneous 
in  such  that  the  intensity  values  are  similar  to  a  pre-determined  probability.  In  the  case  of  roads 
this  probability  is  typically  Gaussian  in  nature  since  the  closer  one  is  to  the  center  the  higher 


the  intensity  and  as  one  moves  towards  the  edge  of  the  road  in  either  direction  the  intensity 
decreases.  The  segment  of  the  road  contour  is  defined  by  Equation  2: 


dv 

dt 


^ik{v)  +  -  log^  + 


V  (T? 


Where  n  is  the  curvature  of  the  contour.  Each  point  along  the  contour  is  encompassed  by  a 


circle  (region)  and  the  value  of  the  mean  and  standard  deviation  of  the  image  pixels  within  this 


circular  region  are  computed,  1  and  S  respectively.  The  mean  and  standard  deviation  of  the  road 


segment  are  defined  as  fj,  and  a  respectively,  and  the  mean  and  standard  deviation  for  the  single 
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point,  P,  which  lies  on  the  eontour,  are  defined  as  fij  and  aj  respeetively.  The  unit  normal  veetor 
to  V  is  n{v).  Onee  the  statistieal  parameters  of  the  point,  P,  are  determined  the  likelihood  ratio 
is  eomputed.  This  determines  whether  the  point  eorrelates  into  the  distribution  for  the  road  or 
is  elassified  as  an  outlier.  Figure  4  shows  the  referenced  regions. 


Fig.  4:  Region  Loeations  [6]  pp.  1194 

Gamba,  Dell’Aequa,  and  Lisini  utilize  a  priori  information  pertaining  to  the  strueture  of  urban 
road  networks  in  eombination  with  a  Fuzzy  Hough  Transform  (FHT)  in  order  to  extraet  the 
desired  road  network  [7].  There  are  three  assumptions  based  on  a-priori  information: 

1)  Urban  roads  align  along  two  direetions,  whieh  are  not  neeessarily  orthogonal 

2)  Roads  are  typieally  eonneeted,  and  dead-end  roads  in  an  urban  setting  are  unlikely 

3)  Road  interseetions  are  dense 

The  first  step  of  this  algorithm  is  to  identify  the  two  main  direetions  of  the  road  network.  These 
direetions  are  used  in  order  to  pre-proeess  the  image  using  an  adaptive  direetional  filter.  The  use 
of  the  FHT  provides  a  histogram  of  the  most  probable  direetions  of  the  initially  extracted  linear 
features,  some  of  whieh  may  not  neeessarily  be  roads  but  most  of  which  are.  If  in  this  histogram 
of  potential  direetions,  there  is  only  one  peak,  then  that  peak  is  considered  twiee,  however  this 
case  is  highly  unlikely.  The  ideal  ease  is  when  the  histogram  contains  two  or  more  evident 
peaks,  the  two  largest  peaks  are  taken  and  these  are  the  direetions  of  the  urban  road  network 
whieh  are  utilized  in  the  direetional  filter.  Onee  the  image  has  been  proeessed  using  the  adaptive 
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directional  filter,  the  FHT  is  again  applied  to  extract  more  accurate  estimates  of  roads.  Some 
errors  can  arise  when  using  the  FHT,  since  it  extracts  only  linear  segments,  curved  areas  require 
a  smaller  window  of  extraction  and  multiple  extractions  may  be  considered.  Also  portions  of 
a  road  may  fail  to  be  detected  due  to  features  in  the  image.  In  this  case  a  method  based  on 
perceptual  grouping  concepts  is  used  in  the  final  processing  of  the  extracted  roads.  This  process 
uses  logical  intuition  to  refine  the  road  network,  concepts  include,  continuity,  collinearity,  and 
proximity.  This  algorithm  is  governed  by  the  selection  of  six  main  parameters  defined  as  follows: 

1)  rud  -  the  maximum  distance  between  two  near  parallel  segments  (fusion  of  segments, 
proximity  concept) 

2)  nie  -  the  maximum  distance  between  the  extremes  of  a  segment  pair  (connection  of 
segments,  continuity  concept) 

3)  rua  -  the  maximum  angle  tolerance  between  segments  (fusion,  collinearity) 

4)  rrig  -  the  maximum  gap  between  extremes  of  potentially  intersecting  segments  and  the 
forecasted  intersection  point  (intersection  proximity) 

5)  rup  -  the  maximum  displacement  between  the  extracted  path  and  its  best  piecewise  linear 
approximation 

6)  mi  -  the  minimum  length  of  a  line  segment 

Amberg,  Coulon,  Marthon,  and  Spigai  propose  a  method  based  on  dynamic  programming 
and  the  Hough  Transform  to  extract  structures  from  high  resolution  SAR  data  in  an  urban 
environment  [8].  The  first  step  in  the  proposed  algorithm  extracts  straight  lines  from  a  local 
region  based  on  the  Hough  Transform  such  that  these  local  regions  are  darker  in  intensity  when 
compared  to  neighboring  regions.  Then  the  algorithm  extracts  curved  features,  first  slightly 
curved  features  and  then  strongly  curved  ones.  Figure  5  outlines  the  steps  in  the  algorithm. 
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Fig.  5:  Algorithm  Outline  [8]  pp.  1785 

The  data  is  pre-processed  using  filters  and  thresholds.  Then  the  data  is  grouped  into  road  and  non¬ 
road  classes  based  on  the  intensity  of  the  pixel,  darker  corresponding  to  roads.  Then  the  large 
straight  segments  of  roads  are  detected  using  the  Hough  Transform.  The  next  step  is  to  identify 
the  segments  of  road  subject  to  slight  curvature.  This  is  done  using  dynamic  programming, 
which  has  the  capability  of  adapting  to  the  given  data.  Dynamic  programming  searches  and 
adds  pixels  which  belong  to  the  same  curve  based  on  the  minimization  of  the  cost  function.  This 
cost  function  is  determined  by  examining  a  tree  of  segments.  Finally  the  results  are  smoothed 
with  local  linear  approximation  and  gap  filling  with  a  distance  criteria.  The  strongly  curved 
sections  of  road  are  approximated  by  a  hyperbolic  model.  The  underlying  assumption  is  that 
two  linear  segments  are  connected  by  a  strongly  curved  segment.  The  hyperbolic  equation  is 
defined  in  Equation  3: 

a(x  —  +  2b{x  —  p)  {y  —  q)  +  c{y  —  q)^  =  1  (3) 

Where  the  center  is  defined  as  (p,  g)  and  the  parameters,  {a,h,c),  define  the  shape  of  the 
hyperbola.  The  center  of  the  hyperbola  is  estimated  as  the  intersection  point  of  the  two  linear 
segments  it  connects.  As  for  the  determination  of  the  shape  parameters,  this  is  done  using  a 
variation  of  the  Hough  Transform.  An  accumulation  matrix  is  built  for  a  pre-defined  vector 
of  values  of  a,  b,  and  c  using  each  pixel  considered  to  lie  on  the  curve.  The  location  in  this 
accumulation  array  with  the  highest  number  of  votes  corresponds  to  the  best  estimate  of  a,  b, 
and  c  based  on  the  pre-defined  vectors. 
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Hu,  Razdan,  Femiani,  Cui  and  Wonka  propose  a  method  based  on  identifying  a  series  of 
footprints,  iteratively  connecting  these  and  then  removing  falsely  identified  paths  using  a  Bayes 
decision  model  [9] .  They  identify  four  main  concerns  when  attempting  to  extract  road  networks 
from  SAR  images: 

1)  Variety  of  roads  (width,  intensity,  shadows)  in  the  same  image 

2)  Differing  intersection  models  (T-shape,  X-shape,  Y-shape) 

3)  Variable  road  width  due  to  sensor  obstructions  such  as  cars  or  building  shadows 

4)  Image  noise  or  physical  connections  between  roads  and  surroundings  causing  leakage 
The  proposed  algorithm  begins  with  some  assumptions  about  roads  and  SAR  images. 

•  A  homogeneous  region  around  a  pixel  is  approximately  rectangular 

•  Roads  are  approximated  as  long  thin  structures  with  a  bounded  width 

•  In  SAR  images  the  intensity  of  roads  are  either  brighter  or  darker  than  the  surroundings 
however  the  intensity  may  not  be  constant  due  to  obstructions  in  the  sensor’s  path  such  as 
vehicles  and  shadows  cast  by  buildings  or  trees 

The  first  step  is  to  identify  the  so  called  footprints,  which  are  polygon  regions  around  an  identified 
road  pixel.  A  road  pixel,  denoted  by  p,  lies  at  the  center  of  a  spoke  wheel.  This  road  pixel  is 
automatically  identified  by  examining  every  pixel  in  the  image  and  its  neighborhood.  There  are 
64  spokes  evenly  spaced  between  0  and  27r.  Starting  at  p,  the  algorithm  moves  along  each  spoke 
and  observes  the  intensity  changes  in  the  next  set  of  pixels.  As  the  distance  from  p  grows  the 
intensity  change  increases.  When  the  difference  in  the  intensity  rises  above  the  standard  deviation 
along  the  spoke,  the  algorithm  terminates  and  has  reached  an  edge  pixel.  This  procedure  repeats 
on  each  of  the  63  remaining  spokes  until  all  edge  pixels  have  been  found.  These  edge  pixels 
are  then  connected  and  the  road  footprint  for  road  pixel,  p,  has  been  found.  In  Figure  6,  a  total 
of  20  footprints  have  been  identified  from  the  image  on  the  left.  Footprints  0-12  correspond 
to  the  road  while  the  remaining  footprints  correspond  to  buildings. 


Approved  for  Public  Release;  Distribution  Unlimited. 
10 


No  =  0 

No=  1 

No  =  2 

No  =  3 

\> 

No  =  4 

No  =  5 

No  =  6 

No  =  7 

No  =  8 

No  =  9 

No=  10 

No  =  12 

No=  13 

No=  14 

No  =  15 

No=  16 

% 
No=  17 

No=  18 

No=  19 

Fig.  6:  Footprints  [9]  pp.  4147 


The  next  phase  is  the  applieation  of  the  toe-finding  algorithm.  This  algorithm  is  intended  to 
trim  the  footprints  and  remove  falsely  identified  branches.  This  algorithm  is  defined  in  the 
following  steps: 

1)  Shift  the  distance  function,  6{i)  so  that  5(0)  is  less  than  the  average  distance.  Where  i  is 
an  integer  corresponding  to  the  spoke  number,  from  0  to  64 

2)  Locate  the  peaks  such  that  the  local  maximum  distance  is  greater  than  the  average  distance 

3)  Remove  the  small  peaks  with  respect  to  the  maximum,  and  all  peaks  less  than  the  average 
are  removed.  Some  peak,  j,  is  removed  if  the  following  holds: 


<  0.25 


(4) 


4)  Merge  close  peaks,  peaks  that  are  within  45  degrees  of  one  another 

5)  Remove  a  peak  if  the  valley  between  the  next  peak  is  not  deep  enough,  average  value 
between  two  peaks,  ii  and  *2  is: 


*2  ~  L  +  1 


(5) 


Then  if  the  following  criteria  is  met,  the  peak  with  the  larger  distance  measure  is  chosen 
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and  the  smaller  is  removed 


^'-^avg 

(5  (zi)  +  5  (Z2)) 


>  0.8 


(6) 


The  thresholds  of,  0.25,  45  degrees,  and  0.8  were  chosen  through  experimentation.  Figure  7 
shows  three  examples  of  footprints  and  the  toe-finding  algorithm  calculations  of  the  distances. 
In  (d)  the  two  peaks  near  the  end  are  merged  since  the  fourth  criteria  in  the  presented  algorithm 
is  met. 


Fig.  7:  Toe-Finding  Algorithm  [9]  pp.  4148 


The  footprints  are  then  classified  into  five  groups: 

1)  Normal  -  two  toes  and  the  turning  angle  is  less  than  45  degrees 

2)  L-Shaped  -  two  toes  and  the  turning  angle  exceeds  45  degrees 

3)  T  -  Shaped  -  three  toes 

4)  X-Shaped  -  four  toes 

5)  Other  -  more  than  four  toes 

With  the  identified  and  trimmed  footprints  the  final  step  of  the  algorithm  is  to  remove  the  falsely 
identified  footprints,  referring  back  to  Figure  6  these  are  the  footprints  labeled  13  through  19 
which  correspond  mainly  to  buildings.  This  is  done  by  implementing  a  Bayes  decision  rule  based 
on  the  estimated  width  of  the  road  at  a  given  vertex. 

Sklarz,  Novoselsky,  and  Dorfan  introduce  the  concept  of  curve  to  curve  fusion  as  opposed 
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to  the  standard  point  to  point  fusion  [10].  A  curve  is  considered  to  be  a  unified  sequence  of 
points.  The  roadmap  is  a  data  structure  consisting  of  nodes  and  road  segments.  These  nodes  are 
the  intersections  of  one  or  more  road  segments.  The  process  of  accepting  a  track  is  as  follows: 

1)  Track  is  accepted  to  the  roadmap  network 

2)  Determine  whether  or  not  the  track  overlaps  with  another  existing  segment 

a)  If  an  overlap  exists,  segment  the  new  track  such  that  the  end  points  correspond  to 
the  overlapped  track  and  proceed  to  fuse 

b)  If  no  overlap  exists,  add  additional  nodes  and  edges  to  the  data  structure 
This  process  is  shown  in  Figure  8  where  the  dashed  line  represents  a  new  track. 


3)  Track  2  Lfap  fusion 


4)  TJp  dat  e  d  Ro  adm  ap- 


Fig.  8:  New  Track  Processing  [10]  pp.  1170 

In  order  to  fuse  two  curve  entities  together  a  distance  measure  is  required,  in  this  case  this  is  done 
using  the  inverse  of  the  joint  covariance  matrix  of  the  two  curves  to  be  fused.  The  dissimilarity 
function,  F(a(r)),  is  defined  to  be  the  distance  between  small  segments  of  the  two  curves  at 
the  location  denoted  by,  r  along  the  fused  curve.  This  function  is  the  squared  statistical  distance 
between  the  two  segments  as  shown  in  Equation  7: 

F(a(r))  =  (CMr))  -  CMr))f  {Pi{t{r))  +  (C,(f(r))  -  CMr)))  (7) 

Where  Ci{t)  and  C2{s)  denote  the  two  curves  to  be  fused  with  covariance  matrices  Pi{t)  and 
P2{s)  respectively.  The  objective  function  for  the  optimally  fused  curve,  a*,  is  defined  in  Equation 
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8: 


fT 

a*  =  arg  min  /  (Ci(t(r))  -  C2{s{t)))^  (Pi(t(r))  +  P2{s{t)))~^  (Ci(t(r))  -  C2(s(r)))  dr 

a{T)  Jo 

(8) 

The  optimal  curve,  a*  is  obtained  using  dynamical  programming  since  there  are  exponentially 
many  solutions  to  obtain  an  alignment  curve  between  two  curves  but  only  one  optimal.  For 
each  point  along  the  curve,  the  cost  is  computed  and  stored.  Then  it  is  a  simple  matter 

to  search  this  stored  data  to  determine  the  minimum.  The  cost  at  the  point  along  the  alignment 
curve,  {i,j)  is  obtained  from  Equation  9: 

D{i,j)  =  min  (^D{i  -  l,j)  + 

+  (9) 

D{i  -  1,  J  -  1)  +  j)) 

Where  5rfe_i,fe  =  \/ (^4  -  with  k  and  k-1  composing  a  segment  along 

the  alignment  curve. 

Shackelford  and  Davis  develop  an  algorithm  based  on  the  capabilities  of  the  Hough  Trans¬ 
form  for  extracting  road  networks  from  high  resolution  imagery  [11].  The  Hough  Transform  is 
explained  in  depth  in  Section  4.3.4  and  is  generalized  here.  This  transformation  takes  a  pixel,  or 
coordinate,  in  the  x-y  space  and  transforms  it  into  a  p  and  9  space.  For  a  single  coordinate  the  value 
of  9  varies  over  a  range  of  -90  to  90  degrees.  If  a  group  of  coordinates  is  shown  to  contribute  to  a 
line,  each  of  these  coordinates  will  share  a  point  of  commonality  in  the  Hough  space,  some  (p,  9) 
pair,  which  identifies  the  line.  Algorithms  for  the  Hough  Transform  typically  consist  of  building 
the  Hough  matrix,  an  x  M  matrix,  corresponding  to  the  lengths  of  the  p  and  9  vectors,  which  is 
an  accumulation  matrix  of  the  (p,  9)  pairs,  identifying  the  largest  values  or  peaks  in  the 
accumulation  matrix,  and  finally  the  ability  to  extract  line  segments  and  fill  gaps  between  similar 
segments  with  the  user  specified  inputs.  The  algorithm  identified  in  [11]  is  based  on  the  Hough 
algorithms  and  iterates  based  on  the  line  segment  size,  decreasing  as  the  iterations  proceed. 
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2,2  Overview 


In  this  report  a  method  for  developing  road  networks  and  eharaeterizing  the  uncertainty  in 
these  estimates  is  developed.  It  is  assumed  that  road  networks  can  be  broken  down  into  two  basis 
functions,  linear  and  quadratic  (ellipse).  The  initial  processing  of  the  data  is  done  by  creating  a 
binary  image  and  extracting  possible  line  segments  using  the  Hough  Transform.  However,  the 
Hough  Transform  does  not  provide  a  measure  of  uncertainty,  therefore  the  Total  Least  Squares 
approach  is  implemented  and  the  Cramer  Rao  lower  bounds  is  derived  from  the  maximum 
likelihood  estimate.  The  Total  Least  Squares  solution  allows  for  an  iterative  estimate  which  is 
updated  in  time  as  additional  measurements  become  available.  The  Least  Squares  solution  will 
be  used  as  an  initial  estimate  for  the  Total  Least  Squares  algorithm.  Once  the  sensor  has  stopped 
receiving  data  the  individual  line  segments  can  be  merged,  extended,  and  blended  to  produce  a 
more  complete  road  network. 

In  Section  4.3.4  the  Hough  Transform  is  presented.  The  mathematical  concepts  behind  the 
transform  for  the  line  are  presented  and  a  simple  example  is  shown  in  order  to  understand 
MatLab’s  built  in  functions  for  the  Hough  Transform.  Then  the  Hough  Transform  is  expanded 
in  order  to  identify  circles  and  ellipses  in  images.  The  complications  with  both  expansions  are 
clearly  explained  and  a  few  algorithms  are  presented  for  the  Circle  Hough  transform.  Finally, 
examples  of  how  to  reduce  the  parameter  space  for  the  Ellipse  Hough  transform  are  presented. 

The  derivations  for  the  Least  Squares,  Total  Least  Squares,  and  the  Least  Squares  ellipse  are 
performed  in  Section  4.4.2.  The  Least  Squares  will  be  an  initial  estimate  for  the  Total  Least 
Squares  algorithm.  Common  derivations  of  the  Total  Least  Squares  make  assumptions  about 
the  noise  associated  with  the  measurements  and  the  system  model,  however  the  case  we  are 
interested  in  involves  a  full  fledged  c  ovariance  m  atrix  a  nd  t  he  T  otal  L  east  S  quares  i  s  derived 
based  on  this  assumption.  The  Least  Squares  solution  for  the  ellipse  is  derived  by  Halir  and 
Flusser  [18]  and  is  presented  and  explained  here. 

Section  3.11  contains  the  derivations  involved  in  the  uncertainty  analysis  for  three  cases.  These 
cases  include  the  straight  line  fit  obtained  by  the  Total  Least  Squares  solution,  the  third  order 
polynomial  fit  obtained  using  MatLab’s  polyfit  function  which  solves  the  problem  in  the  Least 
Squares  sense,  and  finally  the  ellipse  obtained  by  the  Least  Squares  algorithm.  Monte  Carlo 
simulations  and  examples  are  performed  in  order  to  show  the  convergence  of  the  derived 
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solutions. 

In  Section  4.0  the  Graphical  User  Interface  is  thoroughly  explained  and  depicted.  Each 
function  is  explained  and  small  examples  of  important  working  features  are  shown  using  the  first 
available  data  set.  The  required  format  of  the  input  data  as  well  as  the  format  of  the  output  data 
is  outlined  in  detail  as  well. 

The  resulting  interface  is  utilized  along  with  two  available  data  sets  and  in  Section  5.0 
the  results  are  shown.  Since  the  data  was  originally  in  Latitude  and  Longitude  coordinates 
it  is  prudent  to  obtain  the  final  results  in  the  same  coordinate  system.  Results  within  the 
GUI  are  shown  along  with  results  of  externally  utilizing  the  output  data  structure.  The  numerical 
results  of  the  Cramer  Rao  lower  bounds  computations  performed  in  Section  3.11  are  presented 
and  explained  for  each  of  the  provided  data  sets. 

In  Section  5.0  we  make  our  concluding  remarks  about  the  work  done  throughout  this 
report  and  outline  some  potential  future  work  dealing  with  the  implementation  as  well  as 
additional  features  which  may  be  incorporated  at  the  discretion  of  the  user. 
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3.0  METHODS 


3.1  Hough  Transform 

The  Hough  Transform  was  first  introduced  in  1962  by  Paul  Hough  and  was  generalized  in 
1971  by  Richard  Duda  and  Peter  Hart  [12]  and  is  utilized  to  determine  geometric  features  in 
binary  images.  The  most  general  transformation  is  used  to  detect  lines  but  the  concept  can  be 
further  developed  to  detect  circles  [  [12],  [13]]  and  ellipses  [  [13],  [14],  [15]]  in  images.  This 
transformation  uses  a  voting  procedure  in  which  the  algorithm  takes  the  binary  pixel’s  coordinates 
and  computes  the  parameters  required  by  the  geometric  feature.  The  parameters  necessary  vary 
for  each  feature,  a  line  requires  two  parameters,  while  the  circle  three  and  the  ellipse  five.  Thus 
the  accumulation  array  increases  in  dimension  and  the  algorithm’s  computational  complexity 
drastically  increases  in  proportion  to  the  number  of  parameters.  The  Hough  Transform  for  the 
line,  circle,  and  ellipse  are  explored  in  the  next  sections. 

3.2  Line  Hough  Transform 

The  Hough  Transform  eliminates  the  complications  which  arise  with  the  identification  of 
vertical  lines  in  an  image.  This  transform  requires  the  image  to  be  binary  in  nature,  where  white 
pixels  correspond  to  ones  and  black  pixels  correspond  to  zeros.  The  derivation  of  the  Hough 
Transform  requires  basic  trigonometric  properties.  Suppose  we  have  a  line  oriented  as  shown  in 
Figure  9  then  by  defining  the  parameters,  p,  and  6  we  can  derive  the  Hough  Transform.  The 
perpendicular  distance  from  the  origin  to  a  line,  is  denoted  by  p.  The  angle  that  this  distance 
vector  makes  with  the  x-axis  is  6. 

We  note  the  definitions  of  the  cosine  and  sine  functions: 

cos6*  =  -  sin6*  =  -  (10) 

P  P 

Algebraic  substitution  of  a  single  cosine  and  sine  term  in  the  Pythagorean  Theorem: 

-  cos6*  +  -  sin6*  =  1  (11) 

P  P 

It  is  now  a  simple  matter  to  rearrange  Equation  1 1  to  obtain  the  conventional  form  of  the  Hough 
Transform  as  given  by  Equation  12. 

p  =  xcosO  +  ysm.6  (12) 
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Fig.  9:  Parameter  Identifieation 


Now  if  Equation  12  is  rearranged  into  a  slope-intereept  form  we  ean  infer  that  when  6  approaehes 
zero  degrees,  corresponding  to  a  vertical  line,  the  same  issue  occurs  such  that  the  slope  tends 
to  infinity.  However,  using  Equation  12  alleviates  this  problem. 

cos6'  r 

y  = - ^ + 


■ft—  .  .  (13) 

sm  y  sm  y 

The  principle  concept  of  the  Hough  Transform  in  the  line  identification  algorithm  can  be  stated 
as  the  following:  if  two  points  are  collinear  then  they  share  a  pair,  (p,  6)  of  commonality  in 
the  Hough  space.  However,  in  order  to  determine  this  common  pair,  a  Hough  matrix  must  be 
constructed,  this  is  done  by  iterating  over  a  9  range  of  -90  to  90  degrees  for  each  white  pixel, 
which  corresponds  to  the  (x,  y)  coordinates,  in  the  image.  Recall  that  9  is  the  angle  the  p  vector 
makes  with  the  x-axis.  We  can  imagine  that  an  arbitrary  number  of  lines  pass  through  each 
coordinate  at  the  discretized  values  of  9  such  as  in  Eigure  10.  The  figure  shows  two  collinear 
points  with  the  same  discretization  of  9.  The  closest  pair  of  (p,  9)  is  shown.  Each  of  the  dashed 
lines  in  Eigure  10  represents  the  perpendicular  distance  from  the  origin  to  the  solid  line  where 
they  terminate  (two  of  the  distance  vectors  are  hidden  from  view  due  to  the  axes).  Since  p  in  the 
Hough  Transform  is  contained  in  a  finite  length  vector,  the  nearest  discrete  value  of  p,  within 
the  vector,  must  be  determined.  This  is  done  using  a  simple  linear  transformation  as  shown  in 
Eigure  11. 

Thus  to  determine  the  appropriate  RIdx  value  Equation  14  is  utilized,  one  is  added  depending 
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Fig.  10:  Coll  inear  Points 


Fig.  11:  Linear  Transformation  of  p 


on  the  first  index  utilized  in  the  corresponding  language  (MatLab  is  utilized  here). 

RIdx  =  round  (  .r  ^  t-.  (MyR  —  MinR)]  +  1  (14) 

\MaxR-  MinR^  J 

For  each  point  in  the  image,  we  must  iterate  over  the  set  9  range  and  determine  the  corresponding 
Rldx  value  for  each  case.  The  corresponding  index  pair,  (p,  9)  in  the  Hough  matrix  is  then 
incremented  and  this  process  continues  until  each  pixel’s  Hough  Transform  has  been  computed. 
For  each  pixel  in  the  image  with  a  coordinate  {x,y),  the  value  of  p  is  plotted  for  the  set  of 
discrete  values  of  9.  This  is  repeated  for  all  pixels  to  generate  an  image  such  as  the  one  shown  in 
Figure  12.  The  common  point  in  the  Hough  space  is  highlighted  in  the  figure  which  corresponds 
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to  the  intersection  of  the  sinusoidal  curves. 


Hoiisli  Traiisfomi 


-80  -60  -40  -20  0  20  40  60  80 
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Fig.  12:  Hough  Transform 

MatLab  has  a  set  of  built  in  functions  which  can  be  used  to  extract  lines  from  the  Hough  matrix. 
These  functions  consist  of: 

•  hough,  creates  the  Hough  matrix,  options: 

-  ThetaResolution  -  specifies  the  discretization  of  the  9  parameter 

-  RhoResolution  -  specifies  the  discretization  of  the  p  parameter 

•  houghpeaks,  determines  the  (p,  9)  pairs  with  the  highest  number  of  votes,  options: 

-  N  -  specifies  the  number  of  peaks  to  identify  (only  peaks  satisfying  Threshold  will  be 
extracted) 

-  Threshold  -  minimum  number  of  collinear  points  to  extract  a  line 

-  NHoodSize  -  the  region  around  the  identified  peak  which  is  eliminated  upon  identifi¬ 
cation 

•  houghlines,  uses  the  peaks  to  determine  the  associated  line  segments,  returns  a  data  structure 
containing  endpoints  of  line  segments,  options: 

-  FillGap  -  distance  between  two  line  segments  which  will  be  filled  in  order  to  bridge 
gaps 

-  MinLength  -  minimum  allowable  length  of  a  line  segment 
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Once  the  endpoints  and  the  (p,  9)  pair  have  been  identified,  the  eorresponding  data  points  ean 
be  assoeiated  with  the  line  estimate;  this  builds  the  data  elusters  for  the  next  stage. 

3,3  Circle  Hough  Transform 

In  the  line  Hough  there  were  two  unknown  parameters  whieh  defined  the  line.  However  with 
the  eirele  there  are  a  total  of  three  parameters,  the  eenter  eoordinates,  (xo,|/o)  and  the  radius, 
R.  The  generie  equation  for  a  eirele  is  given  to  be: 

{x  -  Xof  +  {y-  yof  =  (15) 

Figure  13  shows  a  eirele  with  the  parameters  identified.  Suppose  a  point,  (a;*,  Ui)  lies  on  the  eirele 

A 


Fig.  13:  Cirele  Parameter  Identifieation 

deseribed  by  the  eenter,  (xo,|/o)  and  the  radius,  R,  we  ean  then  determine  the  eorresponding 
values  of  the  eenter  eoordinates  assuming  that  the  radius  of  the  eirele  is  known.  The  angle,  9 
is  the  angle  between  the  arbitrary  x-axis  loeated  at  xq  and  the  point  lying  on  the  edge  of  the 
eirele.  Using  basie  trigonometry  the  equations  for  the  eenter  eoordinates  are  as  follows: 

xo  =  Xi  —  Rcos9  (16) 

yo  =  yi-  Rsm9  (17) 

Sinee  we  assume  that  R  is  known,  the  aeeumulation  matrix  is  built  in  the  same  manner  as  with 
the  line  Hough  Transform,  with  the  variation  that  9  now  ranges  from  0  to  360  degrees.  Therefore, 
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for  each  point  in  the  image,  the  center  coordinates  are  computed  based  on  the  known  value  of 
the  radius  and  the  discretized  values  of  9.  However,  prior  knowledge  of  the  radius  is  unlikely  and 
therefore  the  Circle  Hough  Transform  (CHT)  becomes  more  complex.  It  is  therefore  necessary 
to  establish  a  vector  of  radius  values  (which  may  be  based  on  the  image  size)  which  can  be 
used  to  search  for  the  optimally  fit  circle.  In  this  manner  for  each  Rj  the  center  coordinates 
are  computed  with  the  discretized  values  of  9.  The  construction  of  the  accumulation  matrix  is 
summarized  as  follows: 

•  Iterate  over  each  value  of  Rj  (dimension  3) 

•  Next  iterate  over  each  coordinate  pixel,  {xi,yi) 

•  Finally  iterate  through  the  discretized  9  vector  and  compute  the  center  coordinates,  (xq,  Dq) 
(dimensions  1  and  2) 

The  accumulation  matrix  is  a  three  dimensional  matrix  as  shown  in  Figure  14.  MatLab  does  not 
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Fig.  14:  Accumulation  Matrix 

contain  any  built  in  functions  for  the  CHT  however  Tao  Peng  has  developed  an  algorithm  which 
solves  many  of  the  complications  which  arise  in  images  such  as  overlapping  circles,  incomplete 
circles,  unknown  radii,  multiple  circles,  and  tightly  spaced  circles.  However  we  begin  with  our 
own  simple  code  to  illustrate  the  concept  and  present  results  from  Peng’s  work  later. 

Using  Equations  16  and  17  it  is  possible  to  construct  a  basic  algorithm  to  illustrate  the  CHT. 
In  reality  the  code  can  be  generalized  to  accept  any  initial  input  based  on  the  user’s  design,  it 
need  not  be  a  binary  image,  however  to  keep  with  the  Hough  Transform’s  convention  we  will 
establish  the  input  to  be  a  binary  image.  We  assume  unequal  and  unknown  radii  for  two  circles 
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as  shown  in  Figure  15. 


Fig.  15:  Unequal  Non-concentric  Circles 

To  begin  with,  the  ranges  need  to  be  specified  on  the  radii  as  well  as  the  potential  location  of  the 
center  coordinates.  Vectors  of  the  center  coordinates  are  established  to  create  indexing  for  the 
accumulator  matrix.  With  the  vectors  specified,  the  computed  values  of  the  center  coordinates 
are  then  rounded  to  the  nearest  index  value  in  the  same  manner  as  the  Rldx  value  in  the  linear 
Hough  Transform.  Then  as  outlined  previously,  we  must  iterate  over  the  range  of  radius  values, 
and  then  each  coordinate  pixel,  and  finally  each  value  of  6.  This  iteration  procedure  builds  the  3- 
dimensional  Hough  matrix  which  contains  the  votes  for  the  parameter  triplet,  {Rj,  xo,yo).  Once 
the  Hough  matrix  is  built,  we  can  then  search  this  matrix  for  the  highest  number  of  votes  and 
attribute  that  3-dimensional  location  to  the  parameters  which  define  a  circle.  Since  we  developed 
this  simplified  version  of  the  CHT  algorithm  with  the  assumption  that  the  circles  have  unequal 
radii,  we  simply  employ  the  houghpeaks  function  after  each  2-dimensional  matrix  is  composed, 
which  then  provides  the  peak  associated  with  each  radius,  Rj.  These  peaks  are  then  stored  into 
a  row  vector  with  the  same  length  as  the  radius  vector  and  when  complete,  we  simply  search 
through  this  vector  for  the  two  highest  values  corresponding  to  the  appropriate  radius  values. 
This  algorithm  is  summarized  below: 

1)  Initialize  radius  and  center  coordinate  vectors 

2)  Set  the  number  of  circles  in  the  image  (assumed  known) 

3)  Iterate  over  the  radius  vector 
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4)  For  each  radius,  iterate  over  each  pixel 

5)  For  each  pixel,  iterate  over  the  Q  range,  0  to  360  degrees  and  compute  the  center  coordinates 

6)  Determine  the  nearest  integer  index  value  for  the  computed  coordinates 

7)  For  each  radius,  using  houghpeaks,  extract  the  single  highest  peak,  corresponding  to  the 
center  coordinates  at  each  radius  value 

8)  Search  through  the  candidate  peaks  vector,  from  previous  step,  and  locate  the  highest  N 
counts,  where  N  is  the  number  of  circles.  This  corresponds  to  the  radii. 

9)  Using  the  corresponding  radius  index  value,  return  to  the  candidate  peaks  vector  and  locate 
the  center  coordinates 
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The  results  of  the  simple  CHT  algorithm  applied  to  the  image  shown  in  Figure  15  is  shown  in  Figure 
16. 


Fig.  16:  Simple  CHT  Algorithm  Results 


As  can  be  seen  from  Figure  16  the  center  coordinates  of  the  circles  are  off  slightly.  This  is  due 
to  the  discretization  and  rounding  in  the  indexing  of  the  computed  values.  Table  I  shows  the 
actual  values  (by  using  the  MatLab  data  cursor)  and  the  CHT  values  for  the  radii  and  the  center 
coordinates.  The  first  circle  corresponds  to  the  smaller  of  the  two  and  the  second  is  the  larger. 

TABLE  I:  CHT  Numerical  Comparison 


Circle  1 

R 

Xo 

Vo 

Circle  2 

R 

Xo 

Vo 

Actual 

48 

146 

176 

Actual 

122 

533 

347 

CHT 

48 

146 

175 

CHT 

122 

533 

345 
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With  this  simplified  algorithm  for  the  CHT  there  are  many  drawbacks.  First,  a  minor  issue  is 
that  the  user  must  specify  the  range  for  the  center  coordinates,  this  is  intended  to  save  space  in 
initializing  the  Hough  matrix.  The  next  problem  is  the  prior  knowledge  of  the  number  of  circles 
in  the  image,  and  then  the  ability  to  properly  identify  these  circles.  The  identification  is  subject 
to  several  potential  scenarios  which  were  previously  described  such  as,  overlaps,  incomplete 
circles  in  the  image,  and  tightly  spaced  circles  (the  potential  to  identify  a  circle  between  the  two 
actual  circles  due  to  the  voting). 

Now  we  will  take  a  look  at  Tao  Peng’s  algorithm  and  its  capabilities  at  handling  the  issues 
identified.  This  algorithm  can  be  found  via  the  Math  Works  website  and  has  the  following  inputs: 

•  img  -  a  grey  scale  image  (should  be  altered  before  input) 

•  radrange  -  the  range  of  radius  values  which  the  algorithm  will  search  over 

•  grdthres  -  gradient  magnitude  threshold,  used  to  remove  the  background  image 

•  fltrdLM  R  -  filter  radius  used  in  the  searching  of  maximum  values  in  the  accumulation  array 

•  multirad  -  parameter  which  allows  for  the  detection  of  multiple  concentric  circles 

•  fltrdaccum  -  filter  to  smooth  the  accumulation  array,  similar  to  standard  MatLab  filters 
Also  the  outputs  of  the  algorithm  are: 

•  accum  -  accumulation  array  from  the  CHT 

•  circen  -  centers  of  all  detected  circles  in  the  image 

•  cirrad  -  radius  values  of  all  detected  circles 

•  dbg_LMmask  -  outputs  the  areas  of  interest  in  detecting  the  local  maxima  in  the  accumu¬ 
lation  matrix 

The  algorithm  first  scans  through  the  image  and  removes  all  background  pixels.  Since  this 
algorithm  requires  a  grey  scale  image  as  opposed  to  the  standard  binary,  the  gradient  magnitude 
will  determine  the  regions  of  higher  intensity,  corresponding  to  brighter  pixels  and  these  will  be 
considered  as  the  desired  data  points  for  the  building  of  the  Hough  matrix.  Once  this  separation 
has  been  achieved,  the  filter  is  applied  to  smooth  the  circle  borders  if  needed  and  also  to  thin 
the  amount  of  pixels  corresponding  to  a  circle.  For  example  if  the  circle  is  completely  filled 
in  with  white  pixels,  the  filtering  allows  a  form  of  edge  detection  in  which  higher  weights  are 
attributed  to  the  points  which  lie  on  an  edge  and  the  inner  points  are  reduced.  The  remainder 
of  the  algorithm  progresses  in  the  same  manner  as  the  simplified  version.  The  accumulation 
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array  for  the  centers  is  built  and  then  searched  for  local  maxima  such  that  the  value  exceeds  the 
threshold.  Once  the  centers  have  been  obtained  the  corresponding  data  points  which  produced 
the  votes  for  this  center  are  then  separated  and  the  accumulation  array  is  built  for  the  radius 
values.  The  multiple  radius  detection  is  possible  in  this  step  with  a  proper  thresholding.  This 
radius  accumulation  array  for  each  center  point  (previously  identified)  is  then  searched  for  the 
one  or  more  local  maxima  depending  on  the  multirad  parameter. 

We  now  present  an  example  from  the  supplied  package.  Figure  17  shows  the  grey  scale  image 
which  contains  overlapping  circles  and  concentric  circles. 


Fig.  17:  Tao  Peng  Example  2 

The  algorithm  is  then  applied  and  we  can  see  the  3-dimensional  accumulation  array  and  the 
produced  results  in  Figure  18 
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Raw  Image  with  Circles  Detected  (center  positioirs  and  radii  irrarked)  ^-D  View  of  the  Accumulation  Array 


(a)  CHT  Results  (b)  Accumulation  Array 

Fig.  18:  Results  from  Tao  Peng’s  CHT  Algorithm 

Another  algorithm  has  been  developed  by  David  Young  whieh  is  mueh  simpler  in  its  imple¬ 
mentation  when  eompared  to  Tao  Peng  and  has  similar  results.  In  Young’s  algorithm  he  saves 
eomputation  time  by  ereating  a  1 -dimensional  array  whose  length  is  N^xNyX  Nr  (multiplieation 
not  matrix  dimensions),  where  N  refers  to  the  seareh  array  of  the  x,  y,  and  r  spaees  respeetively. 
The  eomplete  CHT  is  separated  into  three  funetions  defined  below: 

•  circle  Jiough  -  this  funetion  takes  in  the  image,  and  the  radius  range  and  outputs  the  Hough 
matrix  whieh  in  the  end  is  a  3-dimensional  matrix 

•  circle  Jioughpeaks  -  loeates  the  loeal  maxima  in  the  hough  matrix,  requires  inputs  of  the 
Hough  matrix  and  the  radius  range.  Also  optional  inputs  inelude  the  suppression  neighbor¬ 
hood  for  the  eenter  and  the  radius,  as  well  as  the  number  of  maxima  to  extraet 

•  circlepoints  -  ereates  a  template  of  points  for  eaeh  radius  value,  eentered  at  zero  in  order 
to  build  the  1 -dimensional  array  Hough 

This  algorithm  is  able  to  deteet  multiple  eoneentrie  eireles  as  well  as  overlapping  eireles  and 
re-  quires  roughly  the  same  eomputational  time  as  Peng’s  algorithm  however  it  is  mueh  simpler 
in  its  implementation.  We  show  here  the  results  on  the  same  image  from  Figure  17,  the  green 
eireles  are  the  identified  eireles  by  Young’s  CHT  algorithm. 
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Yoluig’s  CHT 


Fig.  19:  Results  from  David  Young’s  CHT  Algorithm 

Guil  and  Zapata  utilize  a  Fast  Hough  Transform  (FaHT)  in  their  deteetion  algorithms 
for  both  the  eirele  and  ellipse  [13].  The  FaHT  essentially  uses  a  eoarse  quantization  of  the 
initial  parameter  spaee  in  order  to  loeally  identify  regions  of  interest  with  a  high  number  of 
votes  in  the  accumulation  array.  After  these  regions  are  identified  a  finer  level  of  quantization 
is  locally  implemented  and  searched  again  for  the  local  maxima  in  the  accumulation  array. 
This  speeds  up  the  algorithm  since  it  no  longer  searches  the  entire  parameter  space  defined 
with  the  fine  level  of  quantization.  The  CHT  algorithm  is  broken  up  into  two  stages.  In  the 
first  stage  the  candidate  centers  are  located  by  observing  the  intersection  of  the  pixel’s 
gradient  vector.  The  points  considered  in  these  candidate  centers  are  then  stored  to  be 
analyzed  for  the  potential  radius  values.  In  the  second  stage  the  radius  is  determined,  the 
distance  from  the  candidate  center  to  each  stored  pixel  is  computed  and  the  corresponding 
index  in  the  accumulation  array  is  incremented. 
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3,4  Ellipse  Hough  Transform 

The  Hough  Transform  can  be  further  expanded  to  detect  elliptical  shapes  in  images.  However 
as  seen  in  the  CHT  the  parameter  space  increased  from  the  traditional  2-dimensional  space 
to  a  3-  dimensional  space  thus  exponentially  increasing  the  computational  requirements.  In  the 
Ellipse  Hough  Transform  (EHT)  the  dimension  now  increases  to  a  5-dimensional  parameter  space 
defined  by,  the  center  coordinates,  {xo^Dq),  the  major  and  minor  axis,  a’  and  b’  respectively, 
and  finally  the  rotational  angle  of  the  ellipse,  0.  The  computational  complexity  of  the  standard 
EHT  algorithm  has  now  increased  exponentially  as  an  accumulator  array  of  dimension  five  must 
be  built  and  searched  for  the  optimal  set  of  parameters.  Therefore  it  is  common  in  practice  to 
reduce  this  space  as  much  as  possible.  Several  authors  have  separated  the  accumulator  array  into 
two  arrays,  the  first  being  2-dimensional  and  containing  the  potential  center  coordinates  and  the 
remaining  3-dimensional  array  contains  the  major  and  minor  axis  as  well  as  the  rotational  angle 
of  the  ellipse. 

Tsuji  and  Matsumoto  [14]  first  eliminate  long  linear  segments  using  the  traditional  Hough 
Transform.  Then  from  the  remaining  line  segments,  they  search  these  candidates  to  find  parallel 
segments.  The  edge  points  of  these  segments  are  then  averaged  and  the  2-dimensional  accumu¬ 
lator  array  is  then  incremented.  Once  all  the  line  segments  which  failed  to  meet  the  minimum 
length  criteria  have  been  processed  in  this  manner,  the  potential  center  locations  are  extracted 
from  the  2-dimensional  array.  The  edge  points  are  then  stored  in  a  matrix,  however  these  edge 
points  need  not  lie  exactly  on  an  ellipse,  therefore  these  candidate  edge  points  are  then  processed 
to  determine  if  they  do  in  fact  lie  within  an  acceptable  distance  to  the  ellipse. 


Q 


Eig.  20:  Candidate  Ellipse  Points  [14]  pp.  778 
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This  is  done  by  considering  the  angle  between  two  points  as  shown  in  Figure  20.  If  this  angle 
is  90°  ±  S  where  S  is  some  tolerance,  then  the  parr  of  points  lies  on  the  ellipse.  Once  all  pairs 
of  candidate  edge  points  have  been  processed  in  this  manner,  a  least  mean  squares  algorithm  is 
applied  to  determine  the  five  parameters.  The  initial  stages  of  the  algorithm  were  intended  only 
to  identify  the  data  points  which  correspond  to  each  ellipse  and  separate  them  into  their  own 
structures. 

Aguado  and  Nixon  [15]  reduce  the  parameter  space  by  utilizing  gradient  direction  and  tangen¬ 
tial  vectors.  Two  tangential  vectors  are  used  to  determine  the  center  of  the  ellipse.  Extrapolating 
these  vectors  to  their  intersection  point,  T,  and  finding  the  midpoint,  M,  between  the  two  points 
which  lie  on  the  ellipse.  Pi  and  P2,  then  by  extrapolating  the  line,  TM  one  can  locate  the  ellipse’s 
center  as  shown  in  Figure  21. 


Fig.  21:  Ellipse  Center  Focation  [15] 


Using  gradient  information,  the  relation  between  the  ellipse  center  and  the  angular  information 
is  given  by: 


^  ^  y-bp 

x”  X  —  ap 

The  orthonormality  property  of  the  axes  gives  us  the  second  necessary  equation: 

tan(0i  —  p)  tan((;/)2  ~  p)  = 

With  the  following  parameters  defined: 

0i  =  tan"^(^^j  02  =  tan"i(^^j  p  =  tan"i(iF)  ^=5  ^ 


(18) 


(19) 
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Where  ’  denotes  the  first  derivative  and  ”  denotes  the  seeond  derivative  of  the  respeetive 
parameter.  Therefore  rather  than  a  5-dimensional  parameter  spaee  the  EHT  algorithm  has  been 
separated  into  two,  2-dimensional  spaees.  The  first  2-dimensional  spaee  builds  the  aeeumulator 
array  for  the  eenter  eoordinates,  xq,  yo  and  the  seeond  2-dimensional  spaee  builds  the  aeeumulator 
for  the  axes  ratios,  K  and  N. 

In  [13]  the  EHT  parameter  spaee  is  simplified  into  two  2-dimensional  aeeumulation  arrays. 
The  first  of  whieh  is  used  to  identify  the  eenter  eoordinates  in  the  same  method  as  in  [15], 
eaeh  pixel’s  gradient  veetor  is  extrapolated  and  the  interseetion  of  a  multitude  of  gradient  veetors 
defines  a  eandidate  eenter  loeation.  Rather  than  searehing  for  the  remaining  three  parameters,  the 
major  and  minor  axes  and  the  orientation  angle,  the  axes  ratio  is  utilized  to  reduee  the  seeondary 
parameter  spaee  to  a  2-dimensional  spaee.  The  EHT  algorithm  presented  in  [13]  is  also  based 
on  the  East  Hough  Transform  in  whieh  the  axes  ratio  and  orientation  angle  are  searehed  on  a 
eoarse  quantization  level  and  then  the  loeal  maxima  are  obtained  and  a  finer  resolution  area  is 


searehed  to  refine  the  maxima.  The  orientation  angle  and  axes  ratio,  0  and  h  respeetively,  must 
satisfy  Equation  20.  Where  the  axes  ratio  is,  h  =  a? /IP". 

(a;-a;o)cos0-f  (|/-|/o)sin0 

h  =  - - ,  .  ,  ,  . - ^ - -  tan(C  -  0)  (20) 

[x  -  Xq)  sm  (f)+[y-  yo)  cos  0 

The  East  EHT  algorithm  in  stage  two  searehes  the  parameter  spaee  for  the  orientation  angle  and 
axes  ratio  with  the  known  eenter  eoordinates,  {xo,yo),  for  a  given  pixel  in  the  image,  {x,y), 
with  the  angle  of  its  gradient  veetor,  (.  The  axes  ratio  ean  be  transformed  using  the  following 
set  of  equations: 


b  =  ^/h-\x  -  Xo)^  -  {y-  yo)‘^ 


a  =  ^/{x-  Xo)‘^  +  h{y  -  yo)'^ 


In  what  follows  we  will  define  a  method  based  on  the  Eeast  Squares  solution  of  fitting  an 
ellipse  to  a  set  of  data  points  [18].  This  solution  is  explored  in  the  next  seetion  and  supplemented 
with  examples.  The  downside  to  this  implementation  is  that  it  assumes  that  all  data  points  passed 
into  the  algorithm  eorrespond  to  a  single  ellipse. 
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3.5  Maximum  Likelihood  Estimators 

In  this  section  we  will  derive  the  solutions  for  the  standard  Least  Squares,  a  general  form  of 
the  Total  Least  Squares  to  allow  for  a  fully  populated  covariance  matrix,  and  finally  the  Least 
Squares  solution  for  the  ellipse  fitting  problem  as  defined  in  [17]  and  [18].  The  Hough  Transform 
in  the  Section  4.3.4  will  be  utilized  to  identify  groups  of  nearly  collinear  points  in  an  image.  A  best 
fit  line  can  be  obtained  using  the  Total  Least  Squares(TLS)  solution  since  there  is  noise  in  both  the 
X  and  y  directions.  The  TLS  solution  which  will  be  derived  does  not  have  a  closed  form  solution 
and  therefore  it  requires  an  initial  estimate  to  solve  for  a  better  approximation  to  the  linear 
coefficients.  This  initial  estimate  can  be  given  by  a  transformed  version  of  the  Hough  coefficients 
(i.e.  transform  from  p  and  9  to  slope,  m  and  intercept,  b),  however  if  6  is  zero  degrees  then  we 
would  obtain  an  infinite  slope,  which  forces  the  TLS  solution  to  diverge.  Therefore,  rather  than  use 
the  transformed  Hough  Transform  coefficients,  we  will  simply  utilize  a  Least  Squares  estimate  as 
the  initial  guess  for  the  TLS  algorithm. 

3.6  Least  Squares 

In  this  section  we  derive  the  Least  Squares  solution.  The  least  squares  fit  for  a  data  set 
minimizes  the  vertical  distance  between  the  curve  fit  and  the  measurements,  thus  assuming  that 
there  is  noise  in  the  observation.  Figure  22  shows  a  group  of  data  points  with  a  line  fit,  where 
the  minimization  of  the  distance  between  the  measurements  and  the  line  in  the  vertical  direction 
is  the  objective  of  the  Least  Squares  method.  In  order  to  minimize  the  distance  we  begin  by 

establishing  the  error  function  which  is  a  simple  summation  of  vertical  distances  over  each  of 

the  measurements,  M.  A  generalized  version  of  the  measurement  equation  is  given  by: 

y  =  Hx  +  z/  (23) 

Where  H  is  the  M  X  N  matrix  of  measurements,  y  is  a  M  x  1  vector  of  observations,  x  is  a 

X  1  vector  composed  of  the  coefficients  of  the  model  fit,  whether  it  be  linear,  quadratic,  cubic, 
or  a  higher  order  fit  and  the  observation  noise  is  given  by  u.  The  generalized  form  gives  the 
following  error  function: 

F;  =  (Hx-y)'^R-'(Hx-y)  (24) 
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Fig.  22:  Least  Squares  Minimization 


Where  R,  M  X  M,  denotes  the  eovarianee  matrix  of  the  observation  noise.  We  then  need  to 
expand  the  error  funetion  so  that  we  ean  take  the  derivative  to  obtain  the  optimal  solution. 

E  =  x^H^R-^Hx  -  2x^H^R-V  +  y^R^  V  (25) 

Next  we  eompute  the  derivative  of  the  expanded  error  funetion  with  respeet  to  x^  and  set  it 
equal  to  zero  to  obtain  the  optimal  solution. 

^  =  2H^R'^Hx  -  2H^R-V  =  0  (26) 

ctx-' 

Therfore,  the  solution  for  the  estimated  eoeffieients  of  the  data  fit  is  given  by  Equation  27. 

X  =  (H^R-^H)  H^R- V  (27) 

This  is  the  standard  least  squares  solution  used  in  most  regression  models,  where  x  denotes  the 
veetor  of  the  eoeffieient  estimates. 

3,7  Total  Least  Squares 

In  this  section  we  will  derive  the  Total  Least  Squares  solution,  where  the  difference  is  that  the 
Total  Least  Squares  uses  a  model  which  contains  error  in  both  the  H  and  y  terms  and  therefore 
seeks  to  minimize  the  Mahalanbois  distance  between  the  measurements  and  the  curve  fit.  For 
the  case  of  a  simple  linear  regression,  in  which  the  noise  in  x  and  y  are  equal.  Figure  23  is 
shown,  which  displays  the  measurements  and  the  distance  to  be  minimized  (in  this  special  case 
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Fig.  23:  Total  Least  Squares  Minimization 


the  distance  is  perpendicular). 

The  derivation  we  are  interested  is  presented  in  [16],  for  the  case  of  uncorrelated  non- stationary 
noise.  The  noise  is  allowed  to  vary  in  time,  however,  they  are  independent  in  time.  In  other 
words,  the  measurements  are  independent  with  varying  covariance  matrices.  The  model  is  given 
by: 


y  +  i^y 

(28) 

H  -f  z/// 

(29) 

Where  Vy  and  vh  are  the  noise  in  the  observations  and  measurements  respectively.  The  first 
problem  here  is  to  determine  the  suitable  estimates  for  y  and  H  since  the  known  values  are,  y 
and  H.  The  vector  of  known  quantities  is  established  as,  D  =  [H  y]  and  the  vector  of  estimated 
values,  D  =  [H  y].  Thus  the  probability  of  D  given  D  is  given  to  be: 

1 


p  ( D|D 


(30) 


Where  |R|  denotes  the  determinant  of  the  covariance  matrix.  Since  the  measurements  are  as¬ 
sumed  to  be  independent  the  likelihood  function  can  be  rewritten  as  a  product  of  the  individual 
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measurement  PDFs.  Taking  the  negative  of  the  log-likelihood  function  establishes  the  objective 
function  to  be  maximized,  however  the  objective  function  is  constrained  as  follows: 

1  M  ^  ^ 

arg  min  =  2  ~ 

^  i=l 

Subject  to  dfz  =  0  (32) 

Since  the  measurements  are  independent  of  each  other,  the  objective  function  becomes  a  summa¬ 
tion  due  to  the  properties  of  the  log  function.  The  vector  z  is  a  row  vector  defined  as  [x^,  —  1]^ 
and  dj  and  dj  are  the  ith  row  of  their  respective  matrices,  transposed,  giving  a  iV  -f  1  x  1  vector. 
We  must  now  define  the  covariance  matrix,  R.  For  each  measurement  we  have  stated  that  these 
matrices  are  independent  and  unequal.  For  a  single  measurement,  i,  the  covariance  matrix  is: 

pT  p 

^hiVi  ^y^yi 

This  problem  needs  to  be  reformulated  as  an  unconstrained  optimization  problem.  This  can  be 
done  with  the  use  of  Lagrange  multipliers,  which  requires  the  constraint  function  to  be  multiplied 
by  a  scalar  parameter.  A*.  There  are  M  Lagrange  multipliers,  one  for  each  constraint.  The 
unconstrained  objective  function  is  then  modified  by  adding  the  Lagrange  multiplied  constraints: 

J(x)  =  Aid^z  +  A2d^z  +  ...  +  AMd^z+2X](di-di)  (^d^  -  d*)  (34) 

i=\ 

For  a  single  measurement,  i,  the  expanded  objective  function  yields: 


Our  first  concern  is  to  determine  the  optimal  solution  for,  dj.  In  order  to  so,  we  take  the  derivative 

-T 

with  respect  to  dj  and  solve  for  the  Lagrange  multiplier. 

-  R-^di  +  R-^di  =  0  (36) 

d^i 

To  remove  the  inverse  covariance  matrix,  left  multiply  by,  z  R*: 

z'^RiAjZ  =  z^dj  —  z'^di  (37) 
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Recalling  the  constraint,  z  =  0  and  that  the  Lagrange  multiplier,  Aj,  is  a  scalar  quantity,  it  is 
straightforward  to  solve  for  the  Lagrange  multiplier: 


_ 

*  z^Rz 


(38) 


Since  the  multiplication,  ir  RjZ,  will  always  produce  a  scalar  quantity  we  can  simply  divide  this 
term  through.  Substituting  Equation  38  into  36  to  obtain  the  optimal  solution  for  d*: 


RjZ^djZ 

z^Rz 


(39) 


We  will  denote  djZ  as  e*  since  if  performing  the  vector  multiplication  we  see  that  this  is  in  fact 

~T  ^ 

the  error  in  the  measurement,  x  —  y^.  We  can  perform  the  matrix-vector  multiplication  and 
obtain  separate  equations  for  hj  and  y*. 

r  _r  (R/iiftiX  ~  R/iij/i)  Cj 


z^R,-z 


Vi  =  Vi 


(40) 


(41) 


Equations  40  and  41  yield  the  appropriate  estimates  for  the  measurements.  Therefore  we  can 
substitute  the  optimal  solution  found  in  39  into  the  original  constrained  objective  function. 


M 


J  (x)  =  - 


E  4 


2  =  1 


R/Z 

z^RiZ 


R 


-1 


z'^RiZ 


(42) 


Simplifying  and  performing  the  multiplication,  we  obtain  the  following  form  of  the  objective 
function 


2=1 


^  ej  (zR,z 


z'^RiZ 


1  Af  2 

2  z^Rz 

2=1 


h.x-y, 


M 

-Y _ 

2  Y  -  2Rh,y,x  +  Ry^ 


(43) 


ViVi 


Now  we  desire  the  optimal  solution  of  x  from  the  above  equation,  therefore  we  take  its  derivative 
and  set  it  equal  to  zero  as  follows: 


dJ  (x) 

dx 


(i+i) 


y  h* 


y 

i  =  l  X  'Rhihi'X  —  2R/j.y.X  +  Ry^y^ 


a  X 


Ai) 


(i+i) 


R 


^iVi 


X  R/l,;/l,;X  2,'Rli.y-X  T  Ry.y; 


^ViVi 


(44) 


The  denominator  is  denoted  as,  7*  =  x^  R/iih^x  —  2Rhiy^x  +  Ry.y^  and  j  denotes  the  time  step.  We 


must  then  solve  for  the  estimate  at  the  next  time  step,  (j  +  1)  as  follows: 

-1 


M 


Af+1)  -  I 


e?  f  x^-^^ 


R 


hi  hi 


2=1  72  I ^ 


U) 


'jf  ix^^^ 


Vihi 


M 

E 

1=1  7*  I  X 


g2  /  ^(i) 


R 


h'lyi 


.U) 


7?  f  x*--^^ 


(45) 
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No  closed-form  solution  is  available  for  the  updated  estimate,  and  an  initial  estimate  must  be 
given.  This  initial  estimate  can  be  given  by  any  reasonable  estimate  such  as  the  standard  Least 
Squares  solution  or  a  variation  of  the  Total  Least  Squares  with  more  restrictive  assumptions  on 
the  covariance  matrix.  In  the  next  section  we  will  implement  both  the  Least  Squares  and  Total 
Least  Squares  solutions  previously  defined  for  a  simple  line  fitting  problem. 


3,8  LS  and  TLS  Comparison 

First  recall  the  differences  between  the  two  solutions,  the  Least  Squares  minimizes  the  vertical 
distance  whereas  the  Total  Least  Squares  for  the  simulation  which  will  be  presented  here  where 
the  variance  in  x  and  y  are  equal,  minimizes  the  perpendicular  distance.  Equations  46  and  47 
represent  these  errors,  respectively,  which  will  be  computed  for  comparison. 

dy=yi-  (^mxi  +  (46) 

(47) 

+  1 

Equation  47  can  be  obtained  by  projecting  the  vector  from  the  point  of  interest,  i,  to  the  line,  onto 
a  perpendicular  line.  This  is  shown  in  the  next  steps  and  is  also  given  by  WolframMath World 
[20].  Suppose  the  line  of  interest  is  defined  by  the  Equation  48 


ax  +  fey  +  c  =  0 


(48) 


X 

0 

1 

-b 

a  ^  c 

c 

~  b 

a 

b 

Which  can  be  manipulated  into  the  standard  slope-intercept  form  of  a  line.  The  point  of  interest, 
in  which  the  perpendicular  distance  is  to  be  minimized  is  denoted  as,  (xo,|/o)-  Rewriting  in 
slope-intercept  form  and  solving  for  x  and  y  in  terms  of  x  gives  the  vector  equations  in  49. 

(49) 

These  vector  equations  allow  us  to  determine  a  vector  which  is  perpendicular  to  the  line,  v  which 
is  given  by  [a,  hY .  Einally  we  define  the  last  vector,  r,  which  is  a  vector  from  the  point,  (xq,  Dq), 
to  a  point  on  the  line,  (x,  y). 

r 

X  —  Xo 

y-yo 

The  projection  of  r  onto  v  is  a  standard  equation  in  linear  algebra  given  by  Equation  51 

A  _  I''-*’! 

dp 

v 
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(50) 


(51) 


Figure  24  depicts  the  mathematical  equations  in  a  graphical  form. 


ax  +  by  +  c  =  0 


Fig.  24:  Perpendicular  Distance  Point  to  Line  [20] 


Expanding  the  dot  product  in  Equation  51  results  in  Equation  47.  Now  we  present  an  example 
in  which  the  Eeast  Squares  and  Total  Eeast  Squares  solutions  are  used.  A  line  is  defined  by 
its  slope,  m  =  1.243  and  intercept,  b  =  —3.246  with  a  range  of  x  values  of  (—10, 10)  and  a 
step  size  of  0.1.  All  of  the  covariance  matrices,  R,  for  each  measurement  are  assumed  to  be 
equivalent  and  are  defined  as: 


Ris(MxM)  —  0'yi/I(MxM)  —  0.25  X  I(MxM) 


(52) 


2 

^xx 

0 

2 

^xy 

0.25 

0 

0.01 

0 

0 

0 

= 

0 

0 

0 

2 

^yx 

0 

2 

^yy 

0.01 

0 

0.25 

(53) 
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Where  I(mxM)  is  the  M  x  M  identity  matrix  and  M  denotes  the  number  of  measurements.  The  true 
y  values  are  generated  using  the  true  slope,  intereept,  and  x  values.  In  order  to  develop  the 
measurements,  white  additive  Gaussian  noise  is  generated  with  the  above  specified  covariances 
for  each  measurement,  i  =  1,2,  ...,M.  Figure  25  shows  the  true  fit  of  the  line  in  blue,  using 
the  given  slope  and  intercept,  the  generated  measurements  in  black,  the  Least  squares  solution 
in  red,  and  the  Total  Least  Squares  solution  in  green.  The  plot  has  been  zoomed  to  the  region 
shown,  for  a  better  view  of  the  comparison  in  values. 


Least  Squares  vs.  Total  Least  Squares 


Fig.  25:  Least  Squares  and  Total  Least  Squares  Comparison  (zoom) 

Table  II  shows  the  parameter  estimates  for  each  solution  including  the  truth  and  in  addition 
the  root  mean  square  error  is  computed  using  Equations  46  and  47  as  the  errors.  As  shown  in 
the  table  the  error  in  the  estimates  for  the  slope  and  intercept  are  smaller  in  the  Total  Least 
Squares  solution.  Also  the  RMSE  in  the  vertical  distance  is  better  with  the  Least  Squares  while 
the  RMSE  in  the  perpendicular  distance  is  minimal  in  the  Total  Least  Squares,  as  expected. 


TABLE  II:  Parameter  Values  and  RMSE 


Solution 

m 

b 

Error  (m) 

Error  (b) 

RMS  (Vertical) 

RMS  (Perpendicular) 

True 

1.243 

-3.246 

- 

- 

- 

- 

ES 

1.2097 

-3.2683 

0.0333 

0.0223 

1.2014 

0.7805 

TES 

1.2317 

-3.2681 

0.0113 

0.0221 

1.2083 

0.7801 
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3,9  Ellipse  Fitting 

Here  we  present  the  derivation  for  the  Least  Squares  fitting  of  an  ellipse  to  a  set  of  data.  This 
is  presented  by  Fitzgibbon  and  Fiseher  [17]  and  is  further  developed  by  Halir  and  Flusser  [18]. 
This  method  assumes  that  all  of  the  data  passed  into  the  algorithm  is  associated  with  the  ellipse 
and  that  only  one  ellipse  is  present  in  the  data.  Should  there  exist  only  a  segment  of  an  ellipse 
defined  by  the  data,  the  algorithm  will  compute  the  parameters  assuming  an  entire  ellipse  is  to 
be  identified.  We  begin  with  the  basic  equation  for  the  conic  section: 

F{x,  y)  =  ax^  +  bxy  +  cy"^  +  dx  +  ey  +  f  =  Q  (54) 

In  order  to  restrict  the  estimates  of  the  coefficients  to  define  an  ellipse,  the  following  constraint 
must  be  taken  into  account: 

6^  —  4ac  <  0  (55) 

In  terms  of  the  least  squares  sense  we  define  the  vectors: 

n=[x‘^,xy,y‘^,x,y,l]  (56) 

X  =  [a,  b,  c,  d,  e,  ff  (57) 

With  these  vectors  we  can  rewrite  the  original  equation  to,  F(H,  x)  =  Hx.  Now  we  want  to 
minimize  the  sum  of  the  squared  distances  between  the  measurements  and  the  conic  section,  or 
in  this  case  the  ellipse.  Assuming  there  are  M  measurements,  the  constrained  objective  function 
to  minimize  the  distance,  by  determining  the  coefficients  of  the  conic  section,  x,  is  given  by 
Equation  58: 


M 


arg  min 

X 

J  =  F  (h,)^ 

i=l 

(58) 

Subject  to 

4ac  —  b^  =  l 

(59) 

Where  h*  is  the  fth  row  of  the  H  matrix.  The  constraint  has  been  reformulated  into  an  equality 
constraint  by  introducing  an  arbitrary  surplus  variable,  which  is  equated  to  one.  This  is  possible 
because,  for  the  set  of  coefficients  which  define  a  specific  ellipse,  x,  they  can  be  scaled  by  a 
factor,  a,  which  forces  the  equality  constraint,  ax  =  1.  We  prefer  equality  constraints  since  now 
the  constrained  optimization  problem  can  be  solved  using  the  Lagrange  multiplier  method.  In 
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order  to  solve  this  problem  we  must  first  define  a  series  of  matrices,  H  which  is  a  M  x  6  matrix 
of  measurements,  and  a  6  x  6  constraint  matrix,  C. 


H 


0 

0 

2 

0 

0 

0 

x\ 

xm 

yl 

Xi 

yi 

1 

0 

-1 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

yf 

Xi 

yi 

1 

c  = 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

^mVm 

yh 

Xm 

yM 

1  _ 

0 

0 

0 

0 

0 

0 

The  objective  function  is  rewritten  using  the  above  matrices: 


M 


arg  mm 


J  =  ^||Hx| 


i=l 


Subject  to  x^Cx  =  1 


(60) 


(61) 

(62) 


We  multiply  the  constraint  by  the  Lagrange  multiplier,  A  and  append  it  to  the  objective  function 
in  order  to  obtain  the  unconstrained  optimization  function. 


J  =  x^H^Hx  +  A  (1  -  x^Cx) 


(63) 


The  next  step  is  to  differentiate  the  unconstrained  objective  function  with  respect  to  the  coefficient 
vector,  x^: 

^  =  H^Hx  -  ACx  =  0  (64) 

ctx-' 

Solving  for  the  Lagrange  multiplier,  A  by  recalling  the  constraint  equation,  x^Cx  =  1  and  left 
multiplying  the  previous  equation  by  x^  we  get  the  solution  for  the  Lagrange  multiplier: 

A  =  x^Sx  (65) 

Where  the  matrix,  S  =  H^H,  is  termed  the  scatter  matrix.  This  is  a  simple  eigenvalue  problem  in 
which  we  desire  the  minimal  positive  eigenvalue  or  Lagrange  multiplier,  A.  There  exists  up  to  six 
solutions  to  Equation  65.  The  eigenvector  which  corresponds  to  the  minimal  positive  eigenvalue 
is  the  solution  to  the  minimization  problem  and  is  the  best  estimate  of  the  coefficient  vector,  x. 
Halir  and  Flusser  improve  this  algorithm  by  solving  a  set  of  equations  which  are  obtained  by 
partitioning  the  matrices,  S,  H,  and  C  as  well  as  the  coefficient  vector,  x  into  their  linear  and 
quadratic  terms.  They  do  this  because  the  constraint  matrix,  C,  is  singular  and  the  scatter  matrix. 
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S,  can  also  be  singular,  though  very  rarely,  and  only  if  all  of  the  measurements  lie  precisely  on 
the  ellipse.  We  begin  with  the  partitioning  of  the  H  matrix  into  its  linear  and  quadratic  terms  as 
such: 

xf  xiyi  yl  xi  yi  1 

Hi  =  xf  Xiyi  yl:  H2  =  Xi  yi  1  (66) 

.  xIj  XMyn  Um]  I  xm  Vm  ^  _ 


The  remainder  of  the  matrices  are  partitioned  accordingly: 


Solving  Equation  64  using  the  partitioned  matrices  results  in  the  set  of  two  equations: 


SiXi  +  S2X2  —  ACiXi 

(69) 

S^xi  +  S3X2  =  0 

(70) 

We  recall  the  definition  of  the  linear  measurements  matrix,  S3  =  H^H2.  This  matrix  is  singular 
only  when  all  the  measurements  precisely  fit  a  line,  again  this  scenario  for  an  ellipse  is  rare 
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unless  only  points  on  the  long  straight  segments  exist  in  the  measurement  matrix.  Therefore, 
under  the  assumption  that  this  is  not  the  case,  we  can  solve  for  X2  by  taking  the  inverse  of  S3: 

X2  =  -S3  ^S^xi  (71) 

Substituting  Equation  71  into  69  and  since  Ci  is  non-singular  we  get: 

cr'  (Si  -  SaSj^S^)  Xi  =  Axi  (72) 

This  problem  has  been  reduced  from  six-dimensions  to  three,  since  X2  is  now  a  function  of 
Xi.  The  following  set  of  equations  is  easily  solved  by  obtaining  the  eigenvector  and  eigenvalue 
which  pertains  to  the  smallest  positive  number. 


Mxi  =  Axi 


(73) 


xfCiXi  =  1 


X2  =  —S3  ^SoXi 


(74) 

(75) 


Xl 


X  = 


X2 


(76) 


Where  M  =  (Si  —  S2SJ^S^).  This  gives  the  coefficients  of  a  conic  section  constrained  to 
an  ellipse,  however  in  more  practical  scenarios  we  desire  the  parameters  in  terms  of  the  center, 
(xci/o),  major  axis  a',  minor  axis  b'  and  rotation  angle,  0.  The  equation  of  an  ellipse  including 
rotation  is  expressed  as: 

[{x  -  Xo)cos(j)+ {y -yo)sin(j)f  [{x  -  Xq)  sm(j)  +  {y  -  yo)  cos(j)f  _ 

a'2  ^  6'2  -  t  1 
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In  order  to  transform  the  eoeffieients  in  x  into  the  desired  parameters  the  following  equations 
are  used  from  WolframMath World  [19]: 


cd  —  be 

^0  =  7^ - 

—  ac 

ae  —  bd 

w — 
b^  —  ac 

2  (ae^  +  —  2bde  —  acf) 

{h^  —  ac)  \j {a  —  c)^  +  46^  —  (^a  +  c) 

2  (ae^  +  +  fb'^  —  2bde  —  acf) 

—  ac)  —  \l  {a  —  c)^  +  46^  —  (^a  +  c) 

0  for  6  =  0  and  a  <  c. 

for  6  =  0  and  a  >  c. 
^  cot“^  for  6  7^  0  and  a  <  c. 

I  +  I  cot“^  for  6  7^  0  and  a  >  c. 


(78) 

(79) 

(80) 

(81) 

(82) 
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3,10  LS  Ellipse  Example 

In  this  section  we  explore  the  capabilities  of  the  least  squares  ellipse  fitting  algorithm  defined 
by  the  equations  in  the  previous  section.  Since  this  algorithm  is  not  specific  to  the  Hough 
Transform,  it  does  not  require  the  input  to  be  a  binary  image,  instead  it  simply  requires  the  a  set 
of  coordinates  which  are  to  be  considered  part  of  the  ellipse.  The  only  assumption  is  that  all  of 
the  data  points  passed  into  the  algorithm  must  pertain  to  the  ellipse,  as  all  of  these  points  will 
be  considered  when  determining  the  coefficients.  Since  the  algorithm  which  will  be  developed 
in  the  later  sections,  relies  on  the  input  being  a  binary  image,  due  to  the  use  of  the  Hough 
Transform,  we  will  keep  the  same  convention  and  use  an  image  of  an  ellipse  as  our  example 
data.  Figure  26  displays  the  generated  ellipse  (via  the  paint  program)  and  the  results  of  the 
ellipse  Least  Squares  algorithm.  Keeping  in  mind  that  a  generated  ellipse  using  an  image  is  a 
series  of  square  pixels,  the  algorithm  fits  an  ellipse  to  the  set  of  extracted  data  points. 


Least  Squares  Ellipse  Fit 


Fig.  26:  Least  Squares  Ellipse  Fitting  Results 
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3,11  Uncertainty  Analysis 


In  this  section  we  will  present  the  derivation  for  the  Cramer  Rao  lower  bounds  which  provides 
a  measure  of  uncertainty  in  the  coefficients  of  the  linear,  curve,  and  ellipse  fits.  These  derivations 
are  based  on  fully  populated  covariance  matrices.  The  bounds  are  estimated  by  the  inverse  of 
the  Fisher  Information  matrix.  We  begin  by  deriving  the  bounds  for  the  line  fit  and  then  for  the 
blending  functions,  third  order  polynomials  in  x.  We  conclude  the  derivations  with  the  bounds 
for  an  ellipse  defined  by  all  five  parameters.  Each  set  of  derivations  are  accompanied  by  a  Monte 
Carlo  simulation  to  verify  the  convergence  properties  of  the  solutions. 


3.12  Straight  Line 

For  the  straight  line  we  have  a  total  of  three  unknown  parameters,  the  slope,  m,  the  intercept, 
b,  and  the  true  value  of  x,  Xt-  The  measurement  functions  are  given  as: 


y  =  mxt  -\-h  +  Vy 

(83) 

X  =  Xt  +  l'x 

(84) 

The  noise  variables,  Vy  and  are  assumed  to  be  Gaussian  white  noise  and  correlated,  giving 
the  fully  populated  covariance  matrix,  Q. 


Q 


2  2 

^yy 

2  2 

^yx  ^xx 


We  begin  by  determining  the  estimate  of  the  true  value  of  x,  Xt.  This  is  done  using  the  solution 
from  [16]  as  previously  presented.  The  necessary  covariances  are  defined  as  follows: 


R 


hi  hi  — 


R 


hiVi 


^Xiyi  0 


^yryi  —  ^viVi 


(85) 


0  0 

Therefore  the  complete  covariance  matrix,  R,  for  the  Total  Least  Squares  algorithm  is  the 
augmented  matrix  of  the  above  covariances: 

r  ^  2 

'^XiXi  0  '^Xil 

0  0  0 


R, 


R 


hiVi  ^y^yi 


^XiVi  0  ^yrVi 


(86) 


We  now  recall  the  solutions  for  the  estimates  of  hj  and  iji.  Equations  40  and  41  respectively.  Next 
the  equations  are  expanded  by  implementing  the  matrix-vector  multiplications  and  noting  the 
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definitions,  the  covariance  matrix,  R,  the  coefficient  vector,  x  =  [m,  6]^,  and  the  measurements 
and  their  respective  estimated  true  values,  h*  =  [xi,  1],  hj  =  [xj.,  1],  y  =  yt,  and  iji  =  ya-  Thus 
for  a  single  measurement,  i,  we  can  calculate  the  estimated  values  of  x*.  and  yti  as  given  by 
Equations  87  and  88  respectively. 


/^2^2  —  -I- 

^XiXi  ‘^^x^yi  ^  ^ViVi 


Vti  =  Vi 


ma: 


^iVi 


rh?al.^.  —  2mal_...  + 

'^lyt  yiyi 


(87) 


(88) 


Where  we  recall  that,  e*  is  the  error  in  the  measurement  with  the  appropriate  estimates,  e*  = 
mxi  +  b  —  yi.  Note  that  in  Equations  87  and  88,  the  estimate  for  the  slope  and  intercept  must  be 
known  in  order  to  calculate  Xt^  and  ya-  In  our  scenario  these  initial  estimates  are  given  by  the 
basic  Least  Squares  solution.  Since  in  the  Hough  transform,  should  the  value  of  6  be  exactly 
zero  (due  to  the  quantization),  if  we  convert  the  parameters  into  the  necessary  slope  and  intercept 
values,  infinite  values  will  be  obtained.  This  can  be  seen  in  Equations  89  and  90 

cos  9 


m  = 


h  = 


sin^ 

P 


sin  9 


(89) 

(90) 


With  the  initial  estimates  for  the  unknown  parameters,  slope  and  intercept,  one  can  obtain  the 
estimated  true  value,  a;^.,  using  Equation  87.  In  addition  to  the  estimated  true  value  of  x,  the  line’s 
coefficients  are  updated  using  the  Total  Least  Squares  solution.  The  next  step  is  to  determine 
the  uncertainty  in  these  estimates.  The  probability  density  function  for  a  single  measurement, 
{xi,yi),  given  the  estimates  [m,  6,  XtJ  is  given  by: 


p  {xi,yi\m,b,xti )  = 


VL 


Vi 

Xi 


rhxti  +  b 
Xu 


VL 


Vi 

Xi 


rhxti  +  b 
Xu 


(91) 


The  inverse  of  the  covariance  matrix,  is  easily  defined  for  a  2  x  2  matrix: 

.2 


Q 


-1 


cr^  —  cr^ 

^XiXi^yiyi  ^Xiyi 


a, 

—<J. 


2 

MiVi 

2 


-(J. 


^iVi 


a: 


XiXi 


(92) 
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We  expand  the  exponential  term  in  Equation  91  and  write  it  as  a  product  of  two  variables,  ctj 
and  Ki  which  are  defined  in  Equations  93  and  94,  respectively. 


OCi 


K^  = 


2(, 


0-2  0-2  —  0 


iVi  } 


(93) 


^l^x^  i  “  i  +b))  -  (Vi-  ( mxt,  +  b)  )  {xi  -  Xt,)  +  aly^  {xi  -£*,)“ 


(94) 


With  ttj  and  defined  for  each  measurement  we  can  then  rewrite  the  likelihood  function  in  a 
more  contracted  form. 


p{xi,yi\m,b,xti )  = 


MKi 


(95) 


Each  measurement,  i,  is  independent  of  one  another.  Therefore  if  we  assume  there  are  a  total  of 
M  measurements,  the  probability  density  function  for  the  matrix  of  measurements,  [x,  y],  given 
the  parameter  estimates,  [m,b,Xt\,  where  x*  is  now  a  vector  of  estimated  x  values,  is  given  by 
the  product  of  each  measurement’s  probability  density  function: 


M 


P  (x,y|m,6,xj  =  Y[ 


(96) 


=1  \I2ti 


XiXi'^ytyi  ^Xiyi) 

The  Eisher  Information  matrix  is  defined  as  the  negative  expected  value  of  the  Hessian  of  the 
log-likelihood  function  with  respect  to  the  estimated  parameters.  We  define  the  log-likelihood 


function  as,  f  =  In 


p  {x,y\m,b,xt 


,  therefore  the  Eisher  Information  matrix  is  defined  as: 


d^f 

d^f 

d^f 

dmdm 

dmdb 

dmdxt 

d^f 

€l 

d^f 

dbdm 

dbdi) 

dbdxt 

d^f 

d^f 

d^f 

dxtdm 

dxtdb 

dxtdxt 

(97) 


Eirst  the  partial  derivatives  of  /  are  taken  with  respect  to  the  estimated  parameters.  Recall 
that  the  probability  density  function  is  a  product  of  the  individual  measurement’s  and  with  the 
properties  of  the  log  the  product  becomes  a  summation  over  the  M  measurements.  Thus  the 
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partial  derivatives  are: 


M 

Y1  {-XiXt,  +  xl)  +  al^^.  {-2yiXt,  +  2mxl  +  2xtfi^ 


iVl 

tti  -2al.y.  {-Xi  +  Xu)  +  al.^.  {-2yi  +  2mxu  +  2b 


M 

Y1  {~y^  ~  +  ^)  +  ••• 


+  -2|/im  +  2m^xti  +  2m6 


(100) 


Then  the  second  derivatives  which  compose  the  Fisher  Information  matrix  can  easily  be  deter¬ 
mined: 


ni.i)  = 


dmdm 


(101) 


F(1,2)  =  F(2,1)  = 


dmdh 


(102) 


(jp 

F(l,3)  =  F(3,1)  =  =^ai(^-2a2^j,^(-a;i  +  2£iJ  +  a2^,^  [-2y,  +  4mxu  +  2bjj 


(103) 


(104) 


F(2,  3)  =  F(3, 2)  =  ^  a.  +  2<,^m) 

cixtcib 


(105) 


(106) 
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These  equations  then  give  us  an  estimate  of  the  uneertainty  in  the  estimated  parameters.  The 
next  step  is  to  perform  a  Monte  Carlo  simulation  to  show  the  eonvergenee  eharaeteristies  of  this 
estimate.  We  begin  the  simulation  by  ehoosing  a  slope,  intereept,  and  range  of  x  values.  These 
will  be  the  true  simulation  parameters  and  are  speeified  as: 

m  =  -0.4326  b  =  -1.6656  a:  =  -3  :  0.1  ;  3  (107) 

Furthermore  the  eovarianee  matriees,  Q,  for  eaeh  of  the  measurements  is  again  assumed  to  be 
equal  and  is  given  to  be: 

0.5  0.01 

(108) 

0.01  0.5 

Sinee  this  simulation  is  to  determine  the  eonvergenee  eharaeteristies  and  not  the  eapabilities 
of  the  Hough  Transform,  we  will  use  the  Least  Squares  solution  from  Equation  27  where  the 
eovarianee  matrix,  R  is  ayylMxM,  where  M  is  the  number  of  measurements,  whieh  in  this 
simulation  is  61.  Sinee  the  x  truth  was  already  established  the  y  truth  ean  be  ealeulated  using 
the  given  values  for  the  true  slope  and  intereept.  Gaussian  white  noise  is  then  added  to  the 
truth,  whieh  was  speeified  in  Q  and  finally  the  estimated  values  of  x  and  y  can  be  obtained  via 
Equations  87  and  88  respeetively.  A  single  simulation’s  results  are  shown  in  Eigure  27.  Here 
the  truth  is  shown  in  blue,  the  measurements  (truth  with  added  noise)  are  blaek,  and  the  Total 
Eeast  Squares  line  fit  is  given  by  the  red  line. 

The  estimated  values  of  the  slope  and  intereept  from  the  Total  Eeast  Squares  algorithm,  for  this 
single  simulation  run  are: 

rh  =  -0.3775  b  =  -1.5871  (109) 


We  perform  10,000  simulations  to  determine  the  eonvergenee  eharaeteristies  of  the  Eisher  In¬ 
formation  matrix.  The  measure  used  for  eonvergenee  is  the  determinant  of  the  differenee  of 
the  Monte  Carlo  eovarianee  and  the  inverse  of  the  Eisher  Information  matrix.  The  Monte  Carlo 


eovarianee  is  ealeulated  as  a  differenee  of  the  truth  and  the  averaged  estimates.  We  will  denote 
this  eovarianee  as  MCcov  and  is  ealeulated  as: 
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(110) 


Tmth.  Measmeiiients.  and  Estimate 

2 


1 


0  - 

-1- 

-2- 

-3- 

-4- 

-4-3-2-10  1  2  3  4 

X 


■  T  ruth 

Measurements 

■  TLS  Solution 


Fig.  27:  Truth,  Measurements,  Estimate,  TLS  Line  Lit 


Where  the  estimated  values  of  m,  b,  Xi^  are  averaged  after  each  simulation  which  are  denoted 
by,  rfi,  b,  and  Xt^  respectively.  Then  the  convergence  measure  is  given  to  be: 

convergence  =  \MCcov  —  -^”^1  (HI) 

The  Lisher  Information  matrix  is  also  averaged  over  each  simulation.  Ligure  28  shows  the  value 
of  this  convergence  measure  after  each  simulation. 
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Nonnalized  Com’eigeiice 


Fig.  28:  Monte  Carlo  Simulation 


Where  in  Figure  28  we  have  shown  only  a  portion  of  the  total  number  of  simulations  and 
have  normalized  the  eonvergenee  value.  The  estimated  bounds  on  the  parameters,  the  inverse 


of  the  negative  of  the  Fisher  Information  matrix  after  10,000  simulations  is  given  to  be: 


F-^  = 


0.0011 

0 

0.0001 


0 

0.0056 

0.0028 


0.0001 

0.0028 

0.0055 


Next  we  examine  the  estimates  and  their  statistieal  properties.  Eaeh  of  the  estimated  values  of 
m  and  b  are  plotted  in  Figure  29  along  with  the  sigma  ellipses. 

Finally  we  seleet  all  eombinations  of  m  and  b  whieh  fall  within  one  sigma  of  the  average  values 


of  the  respeetive  estimates  and  plot  them,  where  the  range  is  given  as: 


m  = -0.4381  ±0.0336  (112) 

6  = -1.6659  ±0.0719  (113) 
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Slope  and  Intercept  Estimates 
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Fig.  29:  Sigma  Ellipses 

The  lines  defined  by  all  such  coefficients  are  plotted,  in  blue,  along  with  the  true  fit  of  the 
line,  in  red,  in  Figure  30  and  we  can  see  the  uncertainty  in  the  estimates.  As  we  diverge 
from  the  relative  midpoint  of  the  data  range  the  uncertainty  grows.  This  is  as  to  be  expected 
since  the  variance  in  the  y  direction  depends  on  the  variance  of  the  slope,  intercept,  and 
estimated  x  value. 

To  prove  that  the  variance  in  the  y  direction  is  dependent  on  the  variance  of  the  slope,  intercept, 
and  estimated  x  value  we  can  transform  the  Fisher  Information  matrix  from  the  unknown  param¬ 
eters  into  a  variance  in  terms  of  x  and  y.  This  can  be  done  using  the  Jacobian  transformation. 
In  this  transformation  the  Fisher  Information  matrix  is  left  and  right  multiplied  by  the  Jacobian 
of  the  measurement  functions  with  respect  to  the  estimated  parameters.  The  Jacobian  takes  the 
form  of: 


A  = 

dy 

dm 

db 

1 

1 

m 

dx 

dx 

dx 

0 

0 

1 

dm 

db 

dxt 

Therefore  we  perform  the  matrix  multiplication  to  understand  the  growing  variance  in  the  y 
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Fig.  30:  One  Sigma  Slopes/Intercepts 
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This  results  in  the  covariances  taking  the  form  of  Equations  116  through  119 


CT  ^  T  CT 

ViVi  -^ti^mm 


^Xiyi  ^tFmxt  ^bxt  ^^xtxt 


2  _  2 
'^ViXi  ^xyi 


2  _  2 
^XiXi  ^XtXt 


(115) 


(116) 

(117) 

(118) 
(119) 


From  the  above  equations  we  can  see  that  the  variance  in  the  y  direction  depends  on  the  varying 
value  of  Xt^,  so  therefore  as  we  diverge  from  A  =  0  in  either  direction  the  variance  in  y  grows. 


3,13  Algebraic  Fit  Covariance 

In  this  section  we  explain  the  derivation  presented  by  Chernov  and  Lesort  [21]  concerning 
the  covariance  matrix  of  the  weighted  algebraic  fit.  The  algebraic  fit  is  more  commonly  known 
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(120) 


as  the  minimization  of  the  following  summation: 

n 

'^Wi  [P  {Xi,yi-Q)f 

i=l 

In  this  case  x  and  y  are  always  assumed  to  be  the  measurements  and  0  is  the  vector  of  parameters 
which  will  be  defined  for  both  scenarios.  The  weights,  w,  are  assumed  to  be  a  function  of  x,  y, 
and  the  parameter  vector,  0.  The  solution  to  Equation  120  must  satisfy  the  equivalency  to  zero 
of  its  derivative  with  respect  to  the  unknown  parameter  vector,  0.  Using  the  chain  rule  Equation 
121  is  obtained. 

^  +  2^t(;iPiV0P*  =  0  (121) 

The  first  summation  is  discarded  since  we  are  only  solving  with  respect  to  the  leading  order. 
Eurthermore  in  our  case  we  will  assume  that  the  weights,  w,  are  all  equal  to  one.  Therefore  the 
first  term  in  Equation  121  is  nullified  regardless  of  the  leading  order  assumption. 

^WiPiVePi  =  0  (122) 

Since  the  true  points  must  lie  on  the  true  curve,  where  true  is  denoted  by,  x  and  0  ,  using  the 

chain  rule.  Equation  123  can  be  viewed  as  an  additional  constraint  on  the  system. 

(V,P  (x,;  0)  ,  <5x,)  +  (VeP  (x,;  0)  ,  50)  =  0  (123) 

Where  the  angled  brackets  denote  a  vector  product.  We  can  then  substitute  Equation  123  into 
Equation  122  and  solve  for  the  variance  of  the  parameters,  50. 

w,VeP^VeP^SQ  +  w^V^Pf  Sx^VqP,  (124) 

Solving  for  the  variance  term  results  in  Equation  125 

=  -  [5^  wi.VeP^VeP’'^]  [5^  w;,Vxi^^5x, VeP,]  (125) 

The  covariance  can  be  obtained  via  Equation  126. 

Ce  =  E  [5050^]  (126) 
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3,14  3rd  Order  Polynomial 

The  3rd  order  polynomial  is  used  as  the  blending  funetion  between  two  line  segments.  This 
polynomial  takes  the  form  of  Equation  127. 

y  =  Ax^  +  Bx^  +  Cx  +  D  (127) 


We  assume  that  there  is  Gaussian  white  noise  in  both  x  and  y.  Our  desire  is  then  to  determine  the 
uneertainty  of  the  parameter  estimates  determined  in  the  Least  Squares  approach.  The  parameter 
vector  is  defined  as,  0  =  (A,  5,  G,  D).  Using  the  derivations  from  Chernov  and  Lesort  [21],  in 
particular  Equation  125,  the  desired  variance  terms  can  be  obtained.  In  order  to  begin  we  must 
first  define  the  gradient  vectors  with  respect  to  both  the  parameters  and  the  measurements. 

.,3 


VePj  — 


VxPi  = 


dPi 

dA 

dP, 

dB 

dP, 

dC 

dPi 

.  dD  . 

- 

dx 

dy 

x^ 

Xi 


(128) 


“iAxf  +  2Bxi  +  C 
-1 


(129) 


With  the  gradient  vectors  defined,  the  variance  of  the  estimated  parameters  can  easily  be  com¬ 
puted.  Therefore  we  move  onto  a  simulation  to  show  the  results  of  the  derivations  for  the  3rd 
order  polynomial  in  x.  We  begin  by  defining  a  set  of  parameters  and  a  range  on  x. 


yl  =  5  P  =  2  C  =  -l  D  =  A  ;c=-7:0.1:7  (130) 


In  addition  to  the  specified  variables  we  must  define  the  variance  in  the  x  and  y  terms. 

al  =  0.5  al  =  0.75  (131) 

The  plotted  results  of  the  generated  measurements,  in  black,  along  with  the  true  polynomial,  in 
blue,  and  the  estimated  fit,  in  red,  using  the  built  in  MatLab  function,  polyfit,  which  uses  a  Least 
Squares  algorithm  to  determine  the  coefficient  estimates,  are  shown  in  Eigure  31. 

We  then  implement  the  derivations  in  Equation  125  and  obtain  the  following  covariance  matrix 
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3id  Order  Polynomial 


Fig.  31:  3rd  Order  Polynomial  Simulation 
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Figure  32  shows  the  3-sigma  bounded  region  for  the  estimated  polynomial  fit  shown  in  red  in 
Figure  31.  The  plotted  bounds  implemented  here  is  the  upper  and  lower  3-sigma  bounds  on 
eaeh  of  the  parameters  simultaneously.  As  we  defer  from  the  origin  the  bounds  expand  rapidly. 
In  order  to  plot  the  uncertainty  we  consider  two  simple  cases.  These  cases  are  adding  and 
subtracting  the  uncertainty  to  each  of  the  four  estimated  parameters  simultaneously,  these  are 
the  two  cases  shown  in  Figure  32.  The  red  center-line  represents  the  mean  of  the  estimated 
polynomials  and  the  blue  shaded  region  is  the  uncertainty  region  previously  described. 


3-Sigma  Bounded  Region 


Fig.  32:  3rd  Order  Polynomial  Uncertainty 

As  with  the  straight  line  case,  we  also  present  a  Monte  Carlo  simulation.  The  addition  of  the 
noise  along  with  the  polynomial  approximation  and  the  estimation  of  the  bounds  is  performed 
a  total  of  10,000  times.  The  convergence  of  the  solution  is  determined  again  by  a  determinant 
of  the  difference  of  two  matrices. 
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The  first  matrix  is  the  averaged  bounds  eomputed  by  Equation  125  and  the  seeond  matrix  is  the 
averaged  Monte  Carlo  covariance  which  is  calculated  as  in  Equation  133,  however  the  parameters 
now  correspond  to  the  3rd  order  polynomial. 
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(133) 


Where  the  estimated  coefficients  are  averaged  after  each  simulation  which  are  denoted  by,  A, 
B,  C,  and  D  respectively.  Then  the  convergence  measure  is  given  to  be: 

convergence  =  \MCcov  —  (134) 


The  convergence  of  the  normalized  determinant  value  is  shown  in  Eigure  33,  for  only  a  portion 
of  the  total  number  of  simulations. 


Nonnalized  Comergeiice 


Eig.  33:  3rd  Order  Polynomial  Uncertainty  Convergence 

In  addition  to  Eigure  32  we  plot  all  polynomial  fits  with  estimated  parameters  within  3-sigma 
of  the  mean.  The  mean  values  are  given  to  be: 

A  =  3.2252  B  =  1.7460  C  =  41.6646  D  =  7.4202 
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(135) 


From  the  mean  values  we  note  that  there  is  a  large  diserepaney  in  the  value  of  C,  sinee  the 
true  value  was  given  to  be  -1.  The  3-sigma  polynomial  fits  are  shown  in  Figure  34  in  blue 
while  the  polynomial  defined  by  the  mean  values  of  the  parameters  is  shown  in  red. 


3  -  Sigma  Poly  Fits 
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Fig.  34:  3-Sigma  Polynomials 


Also  plotted  in  Figure  34  is  the  true  curve,  in  yellow.  From  this  we  can  see  that  even  though 
there  is  a  large  discrepancy  in  the  value  of  C,  over  the  plotted  range,  there  is  only  a  slight 
difference  in  the  polynomial  fits. 

3.15  Ellipse 

Here  we  derive  the  uncertainty  in  the  parameters  of  the  ellipse  equation  presented  in  Equation 
77  and  is  presented  here  again. 


1  (136) 


We  now  have  five  parameters,  a’,  b’,  xq,  yo  and  0,  in  which  to  characterize  the  uncertainty.  The 
parameter  vector  is  defined  as,  0  =  [a' ,b' ,  Xo,yo,  (j)].  In  order  to  determine  the  uncertainty  we 
need  to  define  the  gradient  vectors  with  respect  to  the  parameters  and  the  measurements.  The 
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gradient  vector  with  respect  to  the  parameters  is  defined  in  Equation  137. 
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(138) 


The  gradient  vector  of  the  ellipse  function  with  respect  to  the  measurements  is  given  by  Equation 
139. 


VxP:  = 


1 

2{(x—XQ)cos(p+{y—yQ)sm.<p)cos(p  ,  2{{x—xo)sm(p—(y—yo)cos(p)sm(p 
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2{{x—xq)  sin  (p—{y—yo)  cos  (p)  cos  <p 

_  dy  _ 

a'3 

62 

(139) 


We  now  want  to  perform  a  simulation  to  show  the  results  of  the  derivations  for  the  ellipse.  The 
following  parameters  are  defined  as  the  true  simulation  parameters. 


TT 


a  =  5  h'  =  2  Xq  =  3  yo  =  2  0=^ 


(140) 


The  variance  in  x  and  y  is  the  same  as  in  the  3rd  order  polynomial  simulation. 


=  0.5  al  =  0.75  (141) 

The  generated  measurements  along  with  the  true  and  estimated  fits  are  shown  in  Eigure  35. 
We  note  here  that  the  algorithm  outlined  in  Section  4.4.2,  derived  by  Halir  and  Elusser  [18]  is 
implemented  in  order  to  determine  the  estimated  parameters. 
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Tmtli.  Measmemeiits.  and  LS  Ellipse  Fit 


X 


Fig.  35:  Ellipse  Simulation 


With  the  generated  measurements  and  the  gradient  vectors,  Equation  125  is  then  utilized  in  order 
to  determine  the  uncertainty  in  the  parameter  vector,  0.  An  example  of  the  covariance  matrix  is 


presented  in  Equation  142. 
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Figure  36  shows  the  uneertainty  in  eaeh  of  the  parameters.  We  have  presented  eaeh  possible 
seenario  for  the  bounds.  For  example,  the  upper  or  lower  bound  on  a  single  parameter,  xq  is 
considered,  and  then  the  filled  region  corresponds  to  the  upper  and  lower  bound  on  the  axis 
parameters,  a’  and  b’.  Therefore  the  upper  and  lower  bounds  on  the  axis  parameters  are  always 
taken  into  account  along  with  each  possible  combination  of  the  remaining  parameters. 


3-Sigiiia  Bounded  Region 


Fig.  36:  Ellipse  Axis  Uncertainty 


We  again  implement  a  Monte  Carlo  simulation  in  which  the  measurements  are  varied,  the 
estimated  ellipse  fit  is  computed,  and  the  covariance  matrix  is  computed  using  Equation  125. 
The  same  method  is  used  to  determine  the  convergence  characteristics  of  the  covariance  solution 
in  which  the  parameters  in  Equation  110  are  replaced  with  the  parameters  in  0  as  shown  in 
Equation  143 
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Where  the  estimated  coefficients  are  averaged  after  each  simulation  which  are  denoted  by,  a! ,  V ,  Xo, 
yo,  and  0  respectively.  Then  the  convergence  measure  is  given  to  be: 

convergence  =  \MCcov  —  -^”^1  (144) 

This  convergence  value  is  normalized  and  the  result  for  a  portion  of  the  simulations  of  the  is 
shown  in  Figure  37. 


Noiinalized  Convergence 


Fig.  37:  Ellipse  Uncertainty  Convergence 

As  with  the  3rd  order  polynomial  fit,  we  also  plot  each  of  the  ellipses  which  is  defined  by 
parameters  within  3-sigma  values  of  their  respective  mean’s.  The  mean  values  of  the  parameters 
are: 

a'  =  4.1178  h'  =  2.3636  Xq  =  2.9998  yo  =  1.9985  0  =  1.0354  (145) 
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Now  each  of  the  ellipses  with  all  five  of  it’s  parameters  lying  within  3 -sigma  values  of  the 
mean  is  plotted  in  Figure  38,  in  blue,  along  with  the  ellipse  defined  by  the  mean  values  of  the 
parameters  in  red. 


3  Sigma  Ellipse  Fits 


Fig.  38:  3-Sigma  Ellipses 
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4,0  ASSUMPTIONS  AND  PROCEDURES 

Throughout  the  course  of  this  section  we  will  explain  the  algorithm  which  has  been  developed 
using  the  theories  from  the  previous  sections.  This  algorithm  automatically  identifies  straight  line 
segments  and  then  through  user  interaction  the  accuracy  of  the  finalized  network  can  be  improved. 
A  graphical  user  interface  (GUI)  has  been  developed  which  allows  the  user  to  easily  interact  and 
manipulate  the  results  at  each  step.  Each  of  the  following  sections  will  explain  the  functions  of 
all  available  buttons  and  as  well  as  how  they  are  implemented. 

4,1  Overview 

The  GUI  consists  of  13  push  buttons,  7  text  boxes,  a  single  slider,  a  radial  button  box  which 
consists  of  2  buttons  (i.e.  only  one  can  be  active  at  a  time)  and  2  plot  areas.  Figure  39  depicts 
the  layout  of  these  features. 
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Fig.  39:  Graphical  User  Interface 

In  what  follows  we  will  display  fragments  of  the  GUI  which  change  after  each  operation  is 
performed  and  the  reader  is  referred  to  Figure  39  as  a  visual  aide  when  needed.  The 
operations  will  be  outlined  as  if  an  actual  database  were  being  proeessed.  For  example  we 
will  start  by  loading  the  database,  then  extracting  the  lines,  manipulating  the  network  through 
the  use  of  the  features  available  within  the  GUI,  and  eventually  extracting  the  network  and 
storing  it  for  later  use.  Figure  40  depicts  a  block  representation  of  the  algorithm  implemented 
during  the  eourse  of  this  section. 


I  Stored  Structure 
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1  Blending 

1  Ellipse  Identification  | 

1  Merging  | 
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Fig.  40:  Algorithm  Overview 
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The  algorithm  begins  by  accepting  a  MatLab  data  structure  which  must  conform  to  a  standard 
set  of  rules. 

1)  Must  be  a  MatLab  data  structure 

2)  Structure  variable  must  begin  with  true,  neglecting  case,  this  field  should  contain  each  of 
the  collected  tracks  (i.e.  trac  ~(1)  is  the  first  track  and  trac  ~(100)  is  the  100th  track) 

3)  Sub-field  containing  the  location  of  the  collected  measurements,  sub-field  contains  loc, 

(trac 

4)  Sub-field  containing  the  covariance  matrices  of  each  of  the  respective  measurements,  sub- 
field  contains  cov,  (trac  .(jf).cov 

The  data  undergoes  some  preliminary  manipulation  which  includes  the  conversion  of  the  co- 
variance  matrix  from  meters  to  degrees  corresponding  to  their  respective  latitudes,  as  well  as 
the  conversion  of  the  measurements  from  latitude  -  longitude  to  pixels  taking  into  account  the 
size  of  the  image  (defaulted  to  1500  x  1500  pixels).  In  addition  to  the  data  conversion  the  data 
from  the  original  structure  is  allocated  into  a  second  structure,  SynData,  which  initializes  the 
following  sub-fields: 

•  SynData.  Location  -  augmented  matrix  of  all  available  measurements  in  both  lat  -  long  and 
pixel  coordinates 

•  SynData.  Covariance  -  augmented  3-dimensional  matrix  of  the  converted  measurement  co- 
variance  matrices 

•  SynData. LastTrack  -  the  number  of  the  last  track  processed 

•  SynData.Xmin  -  minimum  of  the  longitude,  in  meters,  to  offset  the  origin  in  pixel  coordinates 

•  SynData.  Ymin  -  minimum  of  the  latitude,  in  meters,  to  offset  the  origin  in  pixel  coordinates 

•  SynData.Origin  -  the  offset  origin  location  in  latitude  -  longitude 

4.2  Load  Database 

This  push  button  asks  the  user  to  select  an  appropriate  database  to  load.  The  selected  database 
must  conform  to  several  necessary  rules  in  order  to  begin  the  processing.  If  either  of  the  three 
required  fields  is  missing  an  error  box  will  inform  the  user  of  which  field  is  missing.  Once  the 
data  structure  has  been  loaded  there  will  be  a  slight  change  to  the  GUI.  Directly  under  the  Load 
Database  push  button,  there  is  a  string  which  will  state  how  many  tracks  are  available.  In  the 
case  of  the  first  data  set  we  are  working  with,  there  are  4,127  tracks  as  shown  in  Figure  41. 
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Total  Available  Tracks:4,127 


<ri  >l 

400 

Fig.  41:  Load  Database  Change 


In  addition  to  the  notieeable  ehange  to  the  GUI,  the  slider’s  maximum  value  is  updated  sueh 
that  it  eorresponds  to  the  number  of  traeks  available. 

4,3  Line  Extraction 

This  segment  of  the  algorithm  is  the  automated  proeess  whieh  locates  the  lines  in  the  image 
using  the  Hough  Transform,  determines  which  measurements  are  associated  with  these  lines, 
removes  the  measurements  and  covariance  matrices  from  their  initial  fields  in  SynData  and  stores 
them  in  the  respective  line  estimate  field.  Once  the  association  has  been  completed  the  Least 
Squares  and  Total  Least  Squares  solutions  are  computed  and  the  line’s  coefficient  estimates  are 
updated  accordingly.  Several  variables  for  this  process  are  obtained  via  the  GUI.  These  variables 
include;  the  slider  value,  step  size,  distance  threshold,  and  all  of  the  parameters  located  within 
the  Hough  parameters  box  in  the  GUI. 

4.3.1  Slider 

The  algorithm  developed  is  intended  to  be  recursive  such  that  all  tracks  are  not  readily  available 
for  the  given  region.  However,  since  the  test  data  sets  provided  are  composed  of  a  ’’complete” 
set  of  tracks  the  slider  allows  for  a  semi-recursive  view  of  the  data.  A  portion  of  the  data  can 
be  processed  by  setting  the  slider  value  and  then  manipulated  and  stored.  The  stored  data  can 
then  be  loaded  and  the  user  can  continue  working  on  the  data  where  they  left  off.  As  the  slider 
is  moved  the  number  string  below  the  slider  changes  to  denote  the  current  value. 

4.3.2  Step  Size 

This  section  refers  to  the  text  box  to  the  right  of  the  string  of  text.  Step  Size.  The  step  size 
is  a  method  of  segmenting  the  data  in  order  to  process  the  data  recursively.  The  default  value 
of  10  refers  to  10  tracks  processed  at  each  step.  A  step  is  performed  using  the  push  button 
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denoted  by  Step.  This  box  is  one  of  the  neeessary  inputs  prior  to  performing  any  extraction  and 
greatly  impacts  the  Hough  Transforms  extraction  capabilities.  Too  small  a  number  will  result  in 
too  little  data  and  lines  will  not  be  extracted  as  well.  On  the  other  hand,  too  large  of  a  number 
will  result  in  inaccurate  extraction  of  line  segments  from  the  data  (i.e.  the  Hough  will  identify 
undesired  line  segments).  However,  once  the  user  has  established  the  majority  of  the  network 
and  the  Hough  no  longer  identifies  any  segments,  the  user  can  adjust  the  value  to  a  larger  number 
in  order  to  complete  the  data  association  more  rapidly. 

4.3.3  Distance  Threshold 

This  box  allows  the  user  to  specify  the  perpendicular  distance  threshold  used  when  determining 
the  association  of  the  data.  The  input  should  be  given  in  a  pixel  measure,  the  default  being  8 
has  proven  to  be  an  adequate  initial  value.  Once  the  road  network  has  been  further  developed 
the  distance  threshold  can  be  increased  at  the  user’s  discretion  should  there  be  data  which  is  not 
associated  which  should  be. 

4.3.4  Hough  Parameters 

In  this  section  we  will  discuss  the  Hough  Parameters  box  shown  in  Figure  42.  These  are  the 
parameters  which  will  be  passed  into  the  Hough  functions  in  MatLab  in  order  to  extract  the  line 
segments  in  the  image.  These  parameters  were  defined  in  detail  in  Section  and  will  be  reviewed 
here.  Thus  far  the  defaults  for  all  four  of  the  parameters  have  worked  well  and  have  not  been 
altered  during  the  course  of  this  development.  All  of  these  parameters  will  refer  to  a  pixel  metric 
and  not  lat  -  long  or  meters  since  the  Hough  functions  deal  solely  with  an  image. 


|—  Hough  Parameters - 

Rho  Resolution  |  ^  |  Pixels 

Theta  Resolution  ]  i  |  pixels 

Minimum  Length  ]  15  |  pixels 

Fill  Oap  ]  SO  I  Pi«is 


Fig.  42:  Hough  Parameters 
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1 )  Rho  Resolution:  The  rho  resolution  is  the  quantization  of  the  rho  vector,  the  distance  from 
the  origin  to  the  point.  Careful  consideration  must  be  applied  here  as  this  value  along  with  the 
theta  resolution  will  determine  the  ability  of  the  Hough  to  extract  line  segments  appropriately. 

2)  Theta  Resolution:  The  theta  resolution  is  the  quantization  of  the  theta  vector,  the  angle 
of  the  line  created  by  the  point  and  the  x-axis.  The  smaller  the  number  the  larger  the  vector. 
In  later  versions  of  MatLab  the  Hough  function  has  been  changed.  The  alteration  which  has 
been  implemented  deals  with  how  the  theta  option  is  handled.  Rather  than  accepting  a  single 
number  as  the  input  and  having  the  Hough  function  deal  with  the  creation  of  the  vector,  the 
user  is  required  to  input  an  actual  vector  of  values  for  theta.  Therefore  a  version  check  has  been 
incorporated  in  which  the  theta  parameter  issue  is  resolved.  The  theta  vector  remains  between 
-90  and  90  with  a  quantization  given  by  the  user’s  input. 

3)  Minimum  Length:  Recall  that  this  is  the  minimum  allowable  length  of  a  line  segment  to 
be  considered  a  line.  If  the  length,  given  by  the  Euclidean  distance,  is  less  than  the  value  defined 
here  then  the  line  is  not  extracted. 

4)  Fill  Gap:  The  fill  gap  parameter  is  the  distance  between  two  line  segments  with  the  same 
rho  and  theta  value  which  will  be  bridged  if  the  Euclidean  distance  between  the  two  nearest 
endpoints  falls  below  this  threshold. 

4.3.5  Step 

The  implementation  of  this  button  contains  the  extraction  step.  Prior  to  implementing  this 
function  the  user  must  have  changed  each  of  the  previously  mentioned  variables  dealing  with  the 
step  size,  distance  threshold  and  Hough  parameters.  In  addition  to  the  variable  changes  if  desired, 
the  user  must  have  loaded  a  database  to  work  with  which  conforms  to  the  criteria  presented  in 
Section  5.1.  The  first  set  of  tracks  defined  by  the  Step  Size  are  passed  into  the  algorithm  and 
all  of  the  measurements  within  these  tracks  are  manipulated  into  the  new  SynData  structure  and 
converted  as  necessary.  The  conversions  which  occur  include;  converting  original  measurements 
from  latitude  -  longitude  to  meters  and  finally  to  pixels  and  also  the  covariance  matrix  is  converted 
from  the  original  meters  to  degrees  corresponding  to  the  respective  latitude  position.  Equations 
146  to  152  depict  the  necessary  conversion  equations.  The  equations  respectively  represent,  the 
measurement  values  in  meters,  the  measurement  values  in  pixels,  the  respective  covariance  values 
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in  degrees  according  to  their  respective  latitude. 


Xdatai  =  Re  {Latitudeii)  -  Originil)) 

Ydatai  =  Re  {Longitude{i)  —  Origin{2))  *  cos  {Originil)) 

Xindexi  X  -j-  1  -j-  Round 


Xdatai  Xn 


Y, 


index. 


=  2N  —  Round 


datai 


Y 

rr 


_  m 

^XXi  —  ^xx. 


_  m 

^yyi  —  ^yyi 


^xyi  '^■y 


180^ 


‘  {Re  COS  {Latitudeijy 
1802 

{Re  sin  {Latitudei))'^ 
1802 


(146) 

(147) 

(148) 

(149) 
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(152) 


Re  COS  {Latitudci)  sin  {Latitudci) 

In  the  previous  equations  the  following  variables  are  defined: 

•  Re  -  radius  of  the  Earth  at  the  equator  in  meters,  6,371,000  meters 

•  Latitudci  -  latitude  of  the  ith  measurement  in  radians 

•  Longitudci  -  longitude  of  he  emphith  measurement  in  radians 

•  Origin  -  the  latitude  and  longitude,  in  radians,  of  the  first  measurement  used  as  to  shift 
the  data 

•  Xmin  -  the  minimum  value  of  the  Xdata  vector 

•  Y^in  -  the  minimum  value  of  the  Ydata  vector 

•  N  -  size  of  one  of  the  image  cells,  default  value  500 

•  Xinc  -  the  value  of  a  single  pixel  in  the  x-direction,  one  pixel  corresponds  to  20  meters 

•  Yinc  -  the  value  of  a  single  pixel  in  the  y-direction,  one  pixel  corresponds  to  20  meters 

•  a™  , ,  a™ , ,  a^y.  -  respective  covariance  value  in  meters 

In  order  to  allow  the  growth  of  the  image  we  pack  the  initial  image  with  zeros.  For  the  first 
iteration  we  have  a  tiled  image  of  the  following  composition.  A  tiled  image  is  more  desirable 
as  it  significantly  reduces  the  computational  expense  of  the  algorithm. 
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Where  we  have  denoted  the  first  set  of  data  from  the  10  traeks  by  Data.  The  eomplete  image 
size  is  1500  x  1500.  Eaeh  of  the  nine  tiles  is  eheeked  for  data  and  should  there  be  none,  the  tile 
is  ignored  and  the  algorithm  eontinues  to  the  next.  When  data  has  been  found  in  one  of  the  tiles 
the  Hough  Transform  is  then  implemented  and  extraets  up  to  40  lines.  This  number  is  ehosen  to 
be  relatively  high  sinee  the  proeess  will  exit  when  no  line  segment  has  been  found.  The  Hough 
line  extraetion  is  reeursive  in  that  a  single  line  is  identified,  the  measurements  assoeiated  with 
this  line  segment  are  removed  from  the  image  and  then  the  proeess  repeats  until  no  line  segment 
is  identified  by  the  Hough.  The  data  assoeiation  requires  two  thresholds,  one  for  the  distanee 
from  a  point  to  the  line  and  the  seeond  for  the  distanee  from  the  endpoints.  The  first  threshold 
is  obtained  via  the  distanee  threshold  variable  defined  by  the  user,  while  the  seeond  threshold 
is  set  at  5.  This  means  that  the  perpendieular  distanee  from  the  line  to  a  point  must  be  less 
than  8  in  order  to  be  eonsidered.  If  the  first  test  is  passed  then  the  point  is  eheeked  against  the 
endpoints  of  the  line.  Should  the  point  lie  within  the  endpoints  with  the  allowable  toleranee  of 
5,  the  point  will  be  assoeiated  with  the  line  segment  and  removed  from  the  image.  The  endpoint 
toleranee  allows  the  line  to  grow  in  either  direetion  and  therefore  aeross  the  eells  in  the  image. 
Onee  the  Hough  has  finished  extraeting  the  line  segments  and  the  image  has  been  redueed  to 
non-assoeiated  data  the  Total  Least  Squares  algorithm  outlined  in  Seetion  is  implemented.  The 
algorithm  first  eheeks  whether  the  line  is  oriented  in  a  vertieal  or  horizontal  direetion  as  this 
impaets  the  results  of  the  Total  Least  Squares  solution.  A  vertieal  line  is  identified  by  the  Hough 
eoeffieient  of  the  slope  (the  rho  and  theta  ean  be  transformed  to  slope  and  intereept).  Should  the 
absolute  Hough  slope  eoeffieient  be  greater  than  two  the  algorithm  is  reversed  and  solved  as  a 
funetion  of  y  for  the  vertieally  oriented  line.  The  Least  Squares  solution  is  passed  into  the  Total 
Least  Squares  as  the  initial  estimate.  Onee  the  Total  Least  Squares  algorithm  has  finished,  the 
line  segments  as  well  as  the  remaining  non-assoeiated  data  are  plotted  in  the  original  latitude 
-  longitude  eoordinate  system.  This  eompletes  the  proeessing  of  the  10  traeks  and  at  this  time 
the  user  ean  intervene  and  perform  additional  measures  to  modify  the  network  and  inerease  its 
aeeuraey.  These  additional  features  will  be  explained  in  the  later  seetions.  The  results  of  the 
algorithm  are  shown  in  Ligure  43  in  the  first  plot  window  on  the  left,  while  the  seeond  plot 
area  on  the  right  ean  either  display  all  of  the  measurements  eurrently  available  or  the  remaining 
measurements  whieh  were  not  assoeiated  with  the  lines.  In  addition  to  the  obvious  ehanges  in 
the  plot  areas,  the  text  string  above  the  left  plot  has  been  modified  in  order  to  display  the  amount 
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Fig.  43:  10  Tracks  Processed 


of  proeessed  traeks  out  of  the  total  defined  by  the  slider. 

4.3.6  Recursive 

This  button  performs  the  exaet  same  aetion  as  the  Step  button,  however  it  will  proeess  all 
traeks  up  to  the  slider’s  eurrent  value  without  interruption.  It  will  step  through  the  data  as  well, 
for  example  if  there  are  50  traeks  to  be  proeessed  with  a  step  of  10,  the  reeursive  button  will 
proeess  10  traeks  and  show  the  result  of  the  10  traeks  and  then  immediately  eontinue  onto  the 
next  10  and  so  on  until  all  50  are  done  proeessing.  After  the  50  traeks  have  been  proeessed  the 
user  ean  then  manipulate  the  data  as  they  see  fit.  The  reason  for  two  separate  options  is  that 
with  the  Step  button,  the  user  has  more  eontrol  over  the  extraetion  proeess  and  ean  intervene  in 
the  proeessing  of  the  data  to  perform  additional  feature  extraetions.  The  reeursive  funetion  will 
not  allow  any  intervention  during  the  extraetion  proeess  and  the  user  must  wait  until  all  traeks 
have  been  proeessed. 

4,4  Manual  Intervention 

In  this  seetion  the  features  whieh  pertain  to  the  manual  intervention  are  detailed  and  examples 
of  eaeh  are  provided.  The  three  possible  interventions  whieh  ean  be  taken  inelude  merging  of 
two  similar  lines  based  on  the  user’s  diseretion,  the  identifieation  of  an  entire  ellipse,  or  the 
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removal  of  a  previously  defined  feature  whieh  is  undesired  or  not  to  the  user’s  liking  (i.e.  a  line 
identified  by  the  Hough  with  limited  data  may  obseure  the  aeeuraey  of  the  overall  network  and 
ean  be  removed  to  wait  for  additional  data  or  an  adjustment  to  the  Hough  parameters  ean  be 
made). 

4.4.1  Merge 

The  first  of  the  manual  features  allows  the  user  to  seleet  two  similar  lines  whieh  they  believe 
should  be  eombined  to  form  a  single  line.  When  the  button  is  pressed  the  user  will  be  required 
to  seleet  two  lines  by  elieking  near  them  in  the  left  hand  plot  window.  Eaeh  time  the  user  seleets 
a  line,  the  nearest  line  identified  will  turn  from  the  initial  green  eolor  to  red  and  the  user  will  be 
prompted  with  a  eonfirmation  window.  The  eonfirmation  window  eontains  three  options,  ’Yes’, 
’No’,  and  ’Caneel’.  The  ’Yes’  option  will  eontinue  the  proeess  and  the  user  will  be  required  to 
seleet  the  seeond  line  and  eonfirm  the  seleetion  again.  The  response  ’No’  will  revert  the  seleeted 
line  baek  to  green  and  allow  the  user  to  re-seleet  until  they  verify  the  validity  of  the  ehoiee. 
The  third  and  final  option,  ’Caneel’,  will  terminate  the  merge  proeedure  with  no  ehanges  made 
to  the  data  strueture.  The  seleetion  of  two  lines  to  be  merged  is  shown  in  Figure  44,  where  the 
lines  are  identified  in  red  and  the  eonfirmation  box  is  also  shown. 


Fig.  44:  Merge  Seleetion 

The  algorithm  will  then  merge  the  lines  by  eombining  three  sub-fields  eorresponding  to  the 

Approved  for  Public  Release;  Distribution  Unlimited. 

76 


two  line  segments  seleeted.  These  neeessary  sub-fields  include;  Coordinates  (the  measurements 
in  pixel  units),  LatLong  (measurements  in  original  latitude  -  longitude),  and  Covariance  (each 
of  the  respective  covariance  matrices  of  the  measurements  in  degrees).  Once  these  have  been 
combined  the  Total  Least  Squares  algorithm  is  used  to  determine  the  line  fit  of  the  combined 
measurements.  The  resulting  merged  line  of  the  selections  in  Figure  44  is  highlighted  in  Figure 
45. s  The  two  lines  which  were  initially  identified  are  now  removed  from  the  data  structure  and 
replaced  with  a  single  field  corresponding  to  the  merged  line,  s 


Fig.  45:  Merge  Results 


4.4.2  Ellipse 

The  ellipse  button  implements  the  algorithm  defined  by  Halir  and  Flusser  [18].  This  is  similar 
to  the  merge  button  in  that  the  user  is  required  to  define  the  line  segments  which  compose  the 
ellipse.  We  note  here  that  only  a  full  ellipse  is  possible  through  the  use  of  the  algorithm  in 
Section  .  The  user  is  required  to  input  the  number  of  lines  that  compose  the  ellipse  and  then 
identify  each  of  these  lines  in  the  same  manner  as  with  the  merge  algorithm.  After  each  line  is 
identified,  the  confirmation  window  appears,  and  the  user  makes  the  appropriate  choice.  Each 
of  the  lines  is  again  identified  by  changing  the  color  to  red  and  once  all  of  the  lines  have  been 
identified,  the  ellipse  algorithm  takes  all  measurements  in  the  LatLong  sub-field  and  computes  the 
necessary  ellipse  parameters.  In  order  to  compute  the  ellipse  parameters  using  the  measurements 
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in  lat  -  long  coordinates,  it  is  necessary  to  normalize  the  measurements  due  to  the  numerical 
accuracy  of  the  Least  Squares  ellipse  algorithm.  First  the  measurements  are  shifted  such  that 
they  have  a  mean  of  0  in  both  x  and  y,  then  the  shifted  values  are  normalized  by  the  range 
of  their  respective  axis.  These  shifted  normalized  measurements  are  then  passed  into  the  Least 
Squares  ellipse  algorithm.  Figure  46  shows  the  estimated  ellipse.  The  original  lines  will  remain 
in  the  data  structure  however  they  will  be  ignored  for  all  future  operations.  A  sub-field  in  the 
line  structure  labeled,  InEllipse,  will  change  value  from  0  to  1.  This  will  inform  the  algorithm 
to  ignore  any  operations  pertinent  to  the  lines  (i.e.  plotting,  locating,  and  data  association). 


40.686  I - . - 1 - 1 - 1 - rn 

40.684  ■ 

40.682  • 

40.68  ■ 

40.678  -  j 

40.676  ■ 

40.674  - 

40.672  ■ 

-113.92  -113.916  -113.91  -113,905  -113.9 


Fig.  46:  Ellipse  Results 
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The  left  hand  plot  in  Figure  46  shows  the  four  lines  whieh  were  eonsidered  for  the  ellipse 
algorithm.  Once  the  ellipse  has  been  established,  the  measurements  which  were  not  associated 
with  any  line  segments  are  re-checked  to  determine  if  they  should  be  associated  with  the  ellipse. 
This  is  done  by  first  transforming  the  ellipse  to  be  centered  at  the  origin  and  have  a  rotation 
angle  of  zero.  Then  the  remaining  measurements  undergo  the  same  transformation,  in  their 
respective  pixel  units  (easier  to  associate  and  the  lines  used  the  distance  threshold  of  8,  ergo  for 
consistency  the  same  convention  is  used).  This  is  done  in  order  to  simplify  the  optimization 
problem.  In  order  to  transform  the  coordinates  we  refer  to  Equations  154  and  155  from  [24]. 

Xn  =  (x-  Xo)  cos  (t)-{y-  Vo)  sin  0  (154) 

=  -  (a:  -  Xo)  Sin0 -t- (y  -  r/o)  COS0  (155) 

We  denote  the  transformed  points  by  (X„,l^)  and  the  original  measurements  by  {x,y).  The 
center  of  the  ellipse  prior  to  transformation  is  (xo,yo)  and  the  rotation  angle  of  the  ellipse  is 
given  by  0.  With  the  measurements  and  ellipse  transformed  to  the  new  coordinate  system  we 
can  more  easily  compute  the  nearest  point  on  the  ellipse  which  minimizes  the  distance  to  a 
given  measurement.  We  will  be  minimizing  the  Euclidean  distance  between  some  arbitrary  point 
which  lies  on  the  ellipse,  (p,  q)  and  the  measurement  {x,  y).  The  minimization  of  the  Euclidean 
distance  will  yield  the  same  solution  as  the  minimization  of  the  square  of  the  Euclidean  distance. 
Equation  156  represents  the  optimization  problem  to  be  solved,  where  the  solution  is  the  point 
which  lies  on  the  ellipse. 

arg  min  J  =  {p  —  xY  +  {q  —  yY  (156) 

2  2 

Subject  to  ^  ^  =  1  (157) 

In  Equation  156  we  denote  the  semi-axis  values  by  a'  and  b'  respectively.  This  optimization 
problem  can  be  solved  by  the  use  of  the  Eagrange  multiplier  method.  In  this  method  we  append 
the  Eagrange  multiplied  constraint  to  the  cost  function.  The  Eagrange  multiplier  is  given  by,  A. 
The  unconstrained  objective  function  is  then  given  by  Equation  158. 

-  (158) 

In  order  to  solve  Equation  158  we  must  differentiate  with  respect  to  p  and  q,  in  addition  to 
considering  the  constraint.  Performing  the  necessary  calculations  and  simplifying  the  equations 
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arg  min  J  =  {p  -  xf  +  {q  -  yf  +  \  [  ^  ^ 


by  solving  for  p  and  q  we  obtain  the  following  system  of  equations. 


dJ 

(159) 

dp 

P  = 

a'2  +  A 

dJ 

yb'^ 

(160) 

dq 

q  = 

6'2  +  A 

/2  r /2  2  r /2  2/2 

a  0  —  p  0  —  q  a  =0 

(161) 

Substituting  Equations  159  and  160  into  Equation  161  results  in  Equation  162.  In  addition  to 
the  substitutions  we  remove  the  fractions  from  the  equations. 

,,2 


(a'2  +  A)'  (6'2  +  A)' 


+ 


=  0 


(162) 


(a'2  +  A)"  (6'2  -  A)^ 

Multiplying  and  reducing  Equation  162  results  in  a  fourth  order  polynomial  in  A.  The  coefficients 
of  this  polynomial  are  as  follows. 


A": 

1 

(163) 

A': 

2a^  +  2b'^ 

(164) 

A^; 

b'^  +  4a'b'^  +  a'^  -  a'^x^  -  b'Y 

(165) 

A^: 

2a%'^  +  2a%'^  -  2a%'‘^x^  -  2a%’Y 

(166) 

A°: 

a'4^/4  _  ^/2^/4^2  _ 

(167) 

This  solution  will  result  in  four  possible  values  of  A,  two  to  four  of  the  solutions  will  be  imaginary 
since  a  line  can  only  intersect  an  ellipse  at  two  points  (a  single  point  if  it  is  tangent).  Therefore, 
substituting  the  real  solutions  back  into  Equations  159  and  160  will  yield  the  potential  points 
on  the  ellipse.  Then  we  simply  check  the  distances  to  these  points  and  the  minimum  distance 
yields  the  nearest  point. 


4.4.3.  Remove 

This  button  simply  removes  the  selected  feature  and  restores  the  measurements  which  were 
associated  with  the  feature.  Eigure  47  shows  the  removal  of  the  ellipse  which  was  generated 
and  shown  in  Eigure  46.  The  removal  of  the  ellipse  requires  multiple  steps.  This  is  due  to  the 
original  ellipse  being  composed  of  line  segments.  When  the  user  selects  the  ellipse,  the  original 
line  segments  will  be  restored.  Then  should  the  user  wish  to  remove  these  segments  they  may 
do  so  by  again  selecting  the  remove  button  and  individually  selecting  each  line,  as  was  done  in 
Eigure  47. 
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Q  RoadNetworkExtrscfion 

- 

Totel  AvalaUe  Treck$:4,127 


[  Export  Dat^  [  knport  Data  |  Reset  ^  [j  Reinove~j 


I  Display 


Fig.  47:  Removal  of  Ellipse 


4,5  Post  Processing 

The  post  proeessing  box  in  40  refers  to  the  funetions  whieh  should  be  implemented  onee  all 
of  the  available  data  has  been  reeeived  and  proeessed.  These  inelude;  the  use  of  the  blending 
funetion  whieh  allows  the  user  to  seleet  two  lines  and  ereate  a  third  order  polynomial  between 
the  two  lines,  the  trim  and  extend  button  whieh  allows  the  user  to  modify  the  limits  of  a  feature 
with  respeet  to  a  seeond  feature,  the  eomputation  of  the  Cramer  Rao  lower  bounds  estimate 
for  eaeh  of  the  existing  features  (lines,  ellipses,  and  third  order  polynomials),  and  finally  the 
extraet  feature  whieh  allows  the  assoeiation  and  extraetion  of  additional  lines  using  the  remaining 
un-assoeiated  measurements.  Eaeh  of  these  features  is  explored  in  the  next  seetions. 

4.5.1  Blend 

This  button  enables  the  user  to  two  lines,  in  the  same  manner  as  with  the  merge  option,  and 
ereate  a  3rd  order  polynomial,  between  the  seleeted  lines,  in  x  as  defined  by  Equation  168. 

y  =  Ax^  +  Bx^  +  Cx  +  D  (168) 

The  box  next  to  the  blend  button  whieh  is  defaulted  at  0.1  denotes  the  pereentage,  10%,  of  the 
lines  to  eut  and  eonsider  in  the  polynomial  estimation.  Onee  the  user  has  seleeted  and  eonfirmed 
the  two  line  segments  the  algorithm  will  then  determine  the  two  nearest  endpoints  and  trim  the 
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specified  percentage  on  each  line.  The  new  endpoints  will  then  also  correspond  to  the  endpoints 
of  polynomial.  Figure  48  shows  the  polynomial  blending  function  applied  to  two  lines  identified 
in  red.  We  note  here  that  the  polynomial  is  not  updated  when  new  data  becomes  available. 


Fig.  48:  Polynomial  Blend 

Therefore  it  is  ideal  to  utilize  the  blending  after  all  processing  has  been  completed  otherwise  the 
polynomials  will  not  connect  with  the  line  estimates  since  the  line  estimates  are  continuously 
updated. 
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4.5.2  Trim/Extend 


The  Trim/Extend  button  allows  for  a  complete  and  smooth  road  network.  This  function  works 
in  a  particular  order  and  only  certain  cases  have  been  considered.  First  of  all  the  user  must  select 
two  geometric  features,  which  must  be  confirmed  as  in  the  previous  algorithms.  The  first  feature 
selected  will  be  trimmed/extended  with  respect  to  the  second  feature.  The  following  cases  have 
been  considered: 

•  Line  to  Line  -  algorithm  computes  the  intersection  of  the  two  lines  and  then  determines 
the  minimum  distance  of  the  endpoints  of  the  first  line  to  the  intersection  point.  The  point 
which  is  nearer  is  then  replaced  with  the  intersection  point. 

•  Line  to  Curve  -  the  algorithm  computes  the  possible  intersection  points  of  the  line  with  the 
selected  3rd  order  polynomial.  The  distance  between  both  endpoints  and  each  of  the  potential 
intersection  points  is  computed  and  again  the  nearest  intersection/endpoint  combination  is 
identified  and  replaced. 

•  Line  to  Ellipse  -  here  the  algorithm  computes  the  nearest  intersection  point  between  the 
line  and  the  selected  ellipse  and  trims/extends  the  line  to  the  intersection.  If  the  line  already 
intersects  the  ellipse  then  the  endpoint  which  has  passed  through  the  ellipse  is  removed  and 
replaced  with  the  intersection. 

•  Curve  to  Line  -  the  same  as  Line  to  Curve  except  in  this  case  the  Curve  will  be  altered 

•  Curve  to  Curve  -  the  intersection  of  two  3rd  order  polynomials  is  determined  and  trimmed 
accordingly 

Each  of  the  cases  explained  are  presented  in  Figure  49. 
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(a)  Line-Line 


(b)  Line-Curve 


(c)  Line-Ellipse 


(d)  Curve-Line 


(e)  Curve-Curve 


Fig.  49:  Trimmed/Extended  Features 


4.5.3  Uncertainty  Computation 

There  is  no  separate  button  whieh  allows  one  to  eompute  the  Cramer  Rao  lower  bounds 
estimates  of  the  extraeted  features.  Instead,  when  the  ereated  data  strueture  is  saved  using  the 
Export  button,  the  user  is  asked  whether  or  not  they  wish  to  eompute  the  estimate  for  the  Cramer 
Rao  bounds.  Due  to  the  number  of  measurements,  displaying  the  bounds  for  the  ellipse  in  a 
graphieal  manner  is  not  applieable  sinee  as  the  number  of  measurements  grows  the  uneertainty 
in  the  estimate  tends  to  zero,  this  ean  be  seen  in  the  summation  in  the  Fisher  Information  Matrix. 
The  polynomial  blends  whieh  were  ereated  have  a  similar  problem,  where  there  are  either  a  large 
amount  of  measurements  or  very  few  measurements  due  to  the  range  of  the  polynomial.  The 
line  uneertainty  was  eomputed  in  a  series  of  steps  whieh  results  in  the  final  ealeulation  of  the 
uneertainty  in  the  x  and  y  spaee. 

1)  The  data  was  shifted  to  obtain  a  mean  of  zero  in  both  the  latitude  and  longitude. 

2)  Initial  guess  for  the  TFS  algorithm  using  the  previous  slope  and  zero  for  the  intereept. 
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3)  Calculation  of  the  CRLB  using  Equations  101  through  106. 

4)  Jacobian  transformation  to  obtain  the  error  in  the  estimates  of  x  and  y  given  by  Equations 
116  through  119. 

5)  Addition  of  the  measurement’s  covariance  matrix. 

4.5.4  Extract 

This  feature  allows  the  user  to  only  extract  and  associate  data.  The  importance  of  this  feature 
is  such  that  if  the  user  has  already  begun  working  on  blending,  trimming  and  extending  line 
segments,  using  the  Step  button  will  alter  the  lines  back  to  the  original  endpoints  and  the  user 

will  lose  the  work  done.  The  extract  button  allows  only  new  lines  to  be  extracted,  in  the  case 

the  user  has  removed  existing  features  and  wishes  to  modify  the  Hough  Parameters  in  order  to 
detect  smaller  line  segments  or  with  different  orientation.  In  addition  to  extracting  new  segments 
this  can  also  associate  data,  if  the  distance  threshold  is  altered.  Measurements  can  be  associated 
with  a  polynomial  blending  function  based  on  the  nearest  point  on  the  curve  which  is  obtained 
by  solving  the  minimization  of  Equation  169,  where  the  measurement  is  given  to  be  {x^y)  and 
the  point  on  the  polynomial  can  be  defined  as  {xp,  F{xp,  A,  B,  C,  D)). 

Dist  =  \J {x  —  Xp)^  +  {y  —  F{xp,  A,  B,  C,  D))^  (169) 

Differentiation  of  Equation  169  yields: 

dDist  _  -2{x-  Xp)  -  {y-  {Ax^  +  Bxl  +  Cxp  +  D))  {SAx^  +  2Bxp  +  C) 

2\J {x  -  Xp)^  +  {y  -  {Axf,  +  Bx'^  +  Cxp  +  ZJ))^ 

Expanding  Equation  170  and  organizing,  results  in  the  following  equation: 

SA^xl  +  bABxl  +  {AAC  +  2B‘^)  +  {WC  +  ^AD  -  3Ay)  +  ...  (171) 

{C^  +  2CD  -  2By  +  l)  Xp  +  {CD  -  Cy  -  x)  =  0  (172) 

Using  the  roots  command  in  MatEab  gives  the  point  which  lies  on  the  polynomial  and  is  closest 
to  the  measurement  in  question.  We  then  simply  check  the  distance  and  should  it  lie  within  the 
specified  range,  it  is  associated  with  the  polynomial  and  removed  from  the  image. 
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4,6  Stored  Structure 


The  final  box  in  Figure  40  is  the  storing  of  the  processed  data.  Throughout  the  algorithm  a 
MatLab  data  structure  has  been  created  and  continuously  altered  in  order  to  accommodate  the 
necessary  features  and  data  associated  with  such  features.  The  saving  of  this  structure  is  done 
using  the  Export  button  near  the  top  of  the  GUI. 

4.6.1  Export  Data 

This  button  will  performs  two  tasks.  The  first  function  herein  computes  the  Cramer  Rao  lower 
bounds  estimates  for  each  of  the  extracted  features.  The  second  function  is  to  allow  the  user  to 
store  all  the  work  done  during  the  extraction  process  in  a  MatLab  data  structure.  The  structure 
is  outlined  as  follows  with  the  field  names  and  sub-fields  labelled  accordingly. 

•  SynData  -  this  is  the  title  of  the  data  structure. 

-  SynData. Location  -  the  first  field  which  stores  the  measurements  in  a  N  x  4  matrix. 
The  first  two  columns  correspond  to  the  original  latitude  and  longitude  coordinates  and 
the  last  two  columns  correspond  to  the  indexed  y  and  x  values  respectively.  Only  the 
non-associated  measurements  remain  in  this  matrix. 

-  SynData.LastTracks  -  the  number  of  the  last  track  processed. 

-  SynData.Covariance  -  2  x  2  x  N  matrix  of  covariance  matrices  corresponding  to  each 
of  the  non-associated  measurements. 

-  SynData.Ellipse  -  field  of  all  of  the  identified  ellipses,  contains  the  following  sub-fields: 

*  SynData. Ellipse(#). Coordinates  -  contains  the  associated  pixel  measurements. 

*  SynData. Ellipse(#).LatLong  -  latitude  and  longitude  measurements. 

*  SynData.Ellipse(#). Covariance  -  measurement  covariances. 

*  SynData.Ellipse(#).EllipseCoejfs  -  ellipse  coefficients  in  terms  of  parameters,  xq,  yo, 
a’  and  b’. 

*  SynData.Ellipse(#).Phi  -  rotation  angle  of  the  ellipse. 

*  SynData. Ellipse(#).Xtran  -  transformed  x  values  to  be  oriented  with  zero  rotation 
angle  and  centered  at  the  origin  in  terms  of  the  pixel  coordinate  system. 

*  SynData.  Ellipse(#).Ytran  -  transformed  y  values  to  be  oriented  with  zero  rotation 
angle  and  centered  at  the  origin  in  terms  of  the  pixel  coordinate  system. 
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*  SynData.Ellipse(#).CRLB  -  Cramer  Rao  lower  bounds  estimate  matrix  pertaining  to 
each  of  the  parameters  defining  the  uncertainty. 

SynData.PolyBlend  -  field  of  the  blend  functions  with  corresponding  sub-fields: 

*  SynData.PolyBlend(#). Coordinates  -  contains  the  measurements  associated  with  the 
polynomial. 

*  SynData.PolyBlend(#).LatLong  -  measurements  associated  with  the  polynomial  seg¬ 
ment  in  terms  of  latitude  longitude. 

*  SynData.PolyBlend(#).Covariance  -  covariance  matrices  corresponding  to  the  asso¬ 
ciated  measurements. 

*  SynData.PolyBlend(#).Coeffs  -  polynomial  coefficients  from  the  polyfit  function. 

*  SynData.PolyBlend(#).Xspace  -  range  over  which  the  polynomial  exists. 

*  SynData.PolyBlend(#).Fvalues  -  y  values  for  the  polynomial  corresponding  the  xs- 
pace. 

*  SynData.PolyBlend(#).PixelCoejfs  -  coefficients  for  the  polynomial  blend  correspond¬ 
ing  to  the  pixel  space.  Used  in  data  association. 

*  SynData.PolyBlend(#).PixelEnds  -  endpoints  of  the  polynomial  in  the  pixel  space, 
used  in  data  association. 

*  SynData.PolyBlend(#).CRLB  -  the  Cramer  Rao  lower  bounds  estimate  corresponding 
to  the  polynomial  coefficients. 

*  SynData.PolyBlend(#).CRLB_xy  -  the  Cramer  Rao  lower  bounds  estimate  converted 
to  the  x-y  space  using  the  Jacobian  transformation  and  adding  the  original  measure¬ 
ment  covariance. 

SynData.Xmin  -  stores  the  reference  origin’s  x  value  from  the  first  iteration. 

SynData.Ymin  -  stores  the  reference  origin’s  y  value  from  the  first  iteration. 

SynData.Origin  -  the  latitude  longitude  location  of  the  origin,  the  first  measurement 

received  is  considered  the  origin. 

SynData.  Line  Estimate  -  extracted  line  segments  with  sub-fields  listed. 

*  SynData.  LineEstimate(#). Coordinates  -  measurements  associated  with  the  line  seg¬ 
ment  in  terms  of  pixels. 

*  SynData.LineEstimate(#).LatLong  -  measurements  associated  with  the  line  segment 
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in  terms  of  latitude  longitude. 

*  SynData.LineEstimate(#). Covariance  -  covariance  matrices  pertaining  to  the  associ¬ 
ated  measurements. 

*  SynData.LineEstimate(#). Endpoints  -  the  endpoints  determined  by  the  Hough  trans¬ 
formation,  in  pixels. 

*  SynData.LineEstimate(#).Rho  -  Hough  parameter  corresponding  to  the  line. 

*  SynData.LineEstimate(#). Theta  -  Hough  parameter  corresponding  to  the  line. 

*  SynData.LineEstimate(#). Hough  Coefficients  -  coefficients  of  the  line  converted  from 
rho  and  theta  to  slope  and  intercept. 

*  SynData.LineEstimate(#).InEllipse  -  this  tells  several  parts  of  different  algorithms 
whether  the  line  is  found  in  the  ellipse  field. 

*  SynData.LineEstimate(#).LineCoefficients  -  the  estimated  slope  and  intercept  from 
the  TLS  algorithm. 

*  SynData.LineEstimate(#).TrueLat  -  estimated  values  of  latitude,  output  from  the  TLS 
algorithm. 

*  SynData.LineEstimate(#).TrueLong  -  estimated  values  of  longitude,  output  from  the 
TLS  algorithm. 

*  SynData.LineEstimate(#).EstimateEndpoints  -  values  output  from  the  TLS  algorithm 
differ  from  the  Hough  function’s  output. 

*  SynData.LineEstimate(#).E  -  Fisher  Information  Matrix. 

*  SynData.LineEstimate(#).CRLB  -  Cramer  Rao  lower  bounds  matrix  estimate  for  the 
uncertainty  in  the  line  parameters. 

*  SynData.LineEstimate(#).CRLB_xy  -  the  Cramer  Rao  lower  bounds  estimate  con¬ 
verted  to  the  x-y  space  using  the  Jacobian  transformation  and  adding  the  original 
measurement  covariance. 

4,7  Additional  Features 

There  are  a  few  additional  features  in  the  GUI  which  perform  a  few  necessary  tasks  such  as 
resetting  the  GUI  to  the  original  state,  allowing  the  user  to  import  a  pre-existing  data  structure 
and  the  manipulation  of  the  right  hand  plot. 
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4.7.1  Reset 


As  the  button  suggests,  this  will  remove  all  data  from  the  GUI  and  allow  the  user  to  begin 
again  from  serateh.  The  default  values  of  eaeh  of  the  text  boxes  is  also  restored. 

4. 7.2  Import  Data 

This  feature  will  allow  the  user  to  import  an  existing  data  strueture  whieh  was  exported  by 
the  GUI  and  allow  the  user  to  eontinue  working  on  the  extraetion  or  introduee  new  data.  Onee 
the  user  has  imported  the  data  structure  they  can  then  load  a  database  (via  the  Load  Database 
button)  and  continue  the  extraction  process.  The  criteria  of  the  loaded  database  remains  the  same 
as  previously  defined  in  addition  to  the  criteria  that  the  algorithm  will  automatically  prevent  the 
user  from  using  any  tracks  previously  utilized.  For  example  if  the  first  database’s  track  number 
ended  at  500  then  the  second  database’s  first  track  number  must  start  at  501  or  else  the  algorithm 
will  not  allow  the  previous  tracks  to  be  used. 

4.7.3  Image  Options 

The  Image  Options  are  attributed  to  the  second  plot  area  in  the  GUI.  The  All  Measurements 
option  will  display  all  measurements  from  all  tracks  currently  being  processed  while  the  Non- 
Associated  Measurements  will  display  only  the  measurements  which  have  not  been  associated 
with  any  of  the  existing  features.  Both  options  are  shown  in  Figure  50 
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Fig.  50:  Image  Options 
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5.0  RESULTS 


In  this  section  we  import  two  data  structures  available  for  the  purposes  of  testing  the  capa¬ 
bilities  of  the  implemented  algorithm.  The  structure  of  the  supplied  data  is  examined  and  the 
measurements  are  plotted  prior  to  any  transformations.  Each  structure  is  then  loaded  into  the  GUI 
from  Section  4.0  and  the  extraction  process  is  implemented.  The  final  results  of  the 
algorithm  in  graphical  forms  corresponding  to  the  lat  -  long  coordinate  system  are  shown.  In 
addition  to  the  graphical  representation  of  the  network,  examples  of  the  numerical  results 
pertaining  to  the  Cramer  Rao  lower  bounds  estimate  for  a  single  case  of  each  feature  are 

presented. 

5,1  Synthetic  Data 

The  data  which  has  been  supplied  consists  of  simulated  GMTI  tracks  generated  by  the  Air 
Force  Research  Labs  in  Rome,  NY.  Each  track  varies  in  the  number  of  measurements  it  contains, 
however  the  structure  of  each  track  is  consistent.  The  data  structure  is  broken  down  as  follows: 
•  tracks  -  main  field  of  the  structure. 

-  tracks(#).loc  -  Mx  2  matrix  of  latitude  and  longitude  coordinates. 

-  tracks(#).cov  -  4  x  4  xM  covariance  matrices  corresponding  to  each  measurement  of 
the  form: 
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-  tracks(#).vel  -  Mx  2  matrix  of  component  velocities  at  the  same  time  the  measurement 
is  taken. 

-  tracks(#).  update  -  unix  time  representation  of  measurement  time  (seconds  since  January 
1,  1970). 

The  two  supplied  data  structures  follow  the  same  format  with  one  exception.  The  second  data 
set’s  covariance  array  contains  an  extra  matrix  due  to  an  error  in  the  initialization  of  the  array. 
This  extra  matrix  corresponds  always  to  the  first  matrix  and  is  simply  a  matrix  of  zero  values, 
which  is  removed  during  the  pre-processing  stage  of  the  algorithm.  The  first  data  structure 
consists  of  4,127  tracks  and  the  second  contains  1,675,  where  each  track  contains  a  varying 
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amount  of  measurements.  In  Figures  51  and  52  we  show  the  plotted  measurements  from  each 
set  respectively. 
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Fig.  51:  Data  Set  #1  Measurements 
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Fig.  52:  Data  Set  #2  Measurements 


5,2  Data  Set  #1  Results 

In  this  seetion  we  proeess  the  first  available  data  set  whieh  eontains  4,127  traeks  with  a  total  of 
106,936  measurements.  The  GUI  deseribed  in  Seetion  4.0  allows  us  to  alter  the  desired 
variables  pertaining  to  the  line  extraetion  if  need  be.  With  the  first  data  set  the  extraetion 
proeess  took  merely  four  steps: 

1)  An  initial  proeessing  of  approximately  100  traeks  to  provide  a  rough  estimate  of  the  overall 
network  along  with  the  identification  of  the  ellipse. 

2)  Processing  of  remaining  measurements,  4,027  additional  tracks  processed  10  at  a  time  with 
no  user  intervention. 

3)  Merging,  blending,  trimming,  extending,  and  removal  of  appropriate  line  segments. 

4)  Final  association  of  remaining  data  by  increasing  the  distance  threshold  from  8  to  12. 

In  the  next  subsection  we  show  the  graphical  results  in  the  step  by  step  manner  and  then  present 
an  example  of  the  Cramer  Rao  lower  bounds  estimate  for  one  set  of  extracted  features’  parameters 
(i.e.  a  line  segment,  the  ellipse,  and  a  polynomial). 
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5.2.1  Graphical  Results 

The  first  step  was  to  proeess  enough  traeks  until  a  suitable  number  of  lines  were  extraeted  and 
an  overall  estimate  of  the  network  has  been  obtained.  This  took  around  100  traeks  and  Figure 
53  shows  the  rough  estimate  of  the  overall  network.  We  note  that  it  is  not  required  to  obtain  a 
eomplete  representation  of  the  network,  however  it  ean  speed  up  the  proeessing  of  the  traeks. 


Fig.  53:  Data  Set  #1  Phase  One 

This  stage  took  only  a  few  minutes  to  eomplete  and  required  one  ellipse  identifieation  eomposed 
of  four  lines  and  six  line  merging  operations.  Onee  this  was  eomplete  the  remaining  data  was 
proeessed  simultaneously.  The  slider  was  foreed  all  the  way  to  the  right  to  aeeount  for  all 
4,027  remaining  traeks  and  the  step  was  inereased  to  the  number  of  remaining  traeks  in  order 
to  allow  the  algorithm  to  accept  all  measurements  in  the  remaining  tracks.  This  process  took 
approximately  30  minutes  to  complete  using  all  of  the  default  parameters  set  by  the  GUI.  Figure 
54 
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Fig.  54:  Data  Set  #1  Phase  Two 


As  can  be  seen  in  Figure  54,  due  to  the  distance  threshold,  several  lines  were  identified  which 
should  be  merged  with  similar  nearby  line  segments.  These  lines  were  quickly  merged  to  produce 
Figure  55. 


Fig.  55:  Data  Set  #1  Phase  Three 
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The  final  stage  is  to  use  the  blending  and  trim/extend  functions  to  produce  a  fluent  road  network 
which  is  considered  to  be  complete  based  on  the  received  measurements.  The  blending  function’s 
distance  percentage  which  was  trimmed  from  each  of  the  two  lines  varied  from  10%  to  25%  for 
different  scenarios.  Generally  a  higher  percentage  allowed  for  a  longer  less  curved  blend  which 
was  desirable  in  some  cases.  Figure  56  shows  the  final  extraction  with  all  of  the  appropriate 
manual  interactions  implemented. 


Fig.  56:  Data  Set  #1  Final  Extraction  Results 


The  final  results  of  the  extraction  is  29  line  segments,  a  single  ellipse  and  12  polynomial  blending 
functions.  Once  all  of  the  tracks  had  been  processed  it  took  about  15  minutes  to  complete  the 
manual  blending,  trimming,  extending  and  merging  processes  on  the  extracted  features.  Therefore 
to  process  this  entire  data  set  approximately  45  minutes  was  required.  A  total  of  106,915  out  of 
the  106,936  total  measurements  have  been  associated  with  existing  features  in  the  road  network, 
this  corresponds  to  a  99.98%  data  association. 


5.2.2  Uncertainty  Analysis 

Here  we  present  the  results  for  the  uncertainty  pertaining  to  each  of  the  extracted  features  using 
the  equations  derived  in  Section  3.11.  The  Cramer  Rao  lower  bounds  estimates  for  each  of 
the  features  converges  to  very  small  values  due  to  the  summation  across  each  of  the 

measurements  associated  with  the  feature.  For  example,  a  typical  line  contains  anywhere  from 
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measurements  while  the  ellipse  identified  in  this  data  set  eontains  just  over  6,000  measurements 
and  the  polynomials  eontain  on  average  1,500  measurements.  Here  we  will  present  the  results  for 
the  uneertainty  in  eaeh  of  the  estimates  with  numerieal  examples  and  for  the  linear  uneertainty 
we  display  the  Cramer  Rao  lower  bounds  graphieally.  The  linear  uneertainty  ean  be  transformed 
using  the  Jaeobian  transformation  to  obtain  the  uneertainty  in  terms  of  x  and  y  and  supplement 
these  with  their  original  measurement  eovarianee  matriees.  Figure  57  shows  the  Cramer  Rao 
lower  bounds  uneertainty  in  the  estimates  of  the  line’s  parameters  in  red,  whieh  as  ean  be  seen 
these  bounds  are  very  tight. 
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Fig.  57:  Data  Set  #1  Line  CRLB 


We  now  supplement  the  figure  with  a  numerieal  example  of  the  Cramer  Rao  lower  bounds 
estimate  eorresponding  to  eaeh  of  the  parameter  estimates.  This  line  has  4,387  measurements 
assoeiated  with  it. 


(173) 


Using  the  Jaeobian  transformation  we  ean  obtain  the  Cramer  Rao  lower  bounds  in  terms  of  x 
and  y,  this  is  done  using  the  measurements  shifted  sueh  that  they  have  a  mean  of  zero,  sinee 
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this  is  how  the  bounds  were  eomputed.  In  addition  the  Jaeobian  transformed  CRLB  is  added  to 
the  original  eovarianee  matrix  of  the  measurements.  An  example  of  a  measurement  eovarianee 


matrix  is  given  by  Equation  174. 
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The  Jaeobian  transformed  CRLB  is  added  to  the  measurement  eovarianee  defined  in  Equation 


174  to  produee  the  final  eovarianee  in  the  measurement  ineluding  the  eovarianee  in  the  estimate 


given  by  Equation  175. 
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Although  there  is  minimal  ehange  in  the  eovarianee,  the  measurement  eovarianee  is  on  the  order 


of  le-5  while  the  transformed  CRLB  tends  to  be  of  the  order  of  le-9,  sinee  we  are  dealing 
with  latitude  and  longitude  values  a  ehange  of  le-9  is  still  notieeable  in  the  data.  We  now 


present  a  single  line  estimate  from  the  first  data  set  with  several  measurements.  Eigure  58  shows 


the  line  estimate  in  green,  several  measurements  in  blue,  and  eaeh  measurement’s  error  ellipse 
eorresponding  to  one  sigma  value. 
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Eig.  58:  Data  Set  #1  Error  Ellipses  One  Sigma 


Using  the  error  ellipses  we  eompute  the  mean  eovarianee  in  x  and  y.  This  eovarianee  is  then 
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used  to  represent  a  mean  error  in  the  line  using  one  sigma  we  can  see  from  Figure  59  that  this 
encompasses  the  majority  of  the  data  associated  with  each  of  the  line  estimates. 


One— Sigina  Unceilaiiitv'  using  Eiror  Ellipses 


Fig.  59:  Data  Set  #1  Mean  Line  Uncertainty 


The  ellipse  estimate  contains  6,816  measurements  which  causes  the  CRLB  estimate  to  converge 
to  a  very  finite  value,  therefore  we  simply  present  the  numerical  results  here  in  Equation  176. 
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Finally  we  present  the  numerieal  results  for  one  of  the  polynomial  blends  whieh  eontains  1,550 
measurements  given  by  Equation  177. 
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5,3  Data  Set  #2  Results 

In  this  seetion  we  proeess  the  seeond  data  set  whieh  eontains  1,625  traeks  with  a  total  of 
88,685  measurements.  The  seeond  data  set  as  shown  in  Figure  52  was  expeeted  to  be  more 
eomplex  in  all  aspeets  of  the  proeessing.  Not  only  were  more  lines  extraeted  but  additional 
attempts  at  extraeting  were  taken  by  varying  the  parameters  in  the  Hough  seetion  of  the  GUI. 
The  extraetion  proeess  took  a  total  of  seven  steps  outlined  as  follows: 

1)  The  initial  proeessing,  wherein  no  intervention  was  taken,  the  first  400  traeks  were  pro- 
eessed  with  a  step  of  10  and  then  the  remaining  traeks  were  proeessed  as  a  whole. 

2)  The  first  eleaning  stage,  removal  of  several  areas,  merges  of  a  few  lines  and  polynomial 
blends. 

3)  The  seeond  extraetion  stage,  altered  the  parameters  to  attempt  to  extraet  smaller,  finer  line 
segments. 

4)  Seeond  eleaning  stage,  additional  blends  ereated. 

5)  Third  extraetion  stage,  additional  finer,  smaller  segments  identified. 

6)  Final  eleaning  stage. 

7)  One  last  assoeiation  to  assoeiate  neeessary  data. 

In  the  next  subseetion  we  show  the  graphieal  results  in  the  step  by  step  manner  and  then  present 
an  example  of  the  Cramer  Rao  lower  bounds  estimate  for  one  set  of  extraeted  features’  parameters 
(i.e.  a  line  segment,  the  ellipse,  and  a  polynomial). 
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5.3.1  Graphical  Results 

The  second  data  set  was  completely  processed  without  manual  intervention  (i.e.  merging  or 
ellipse  finding).  This  data  set  as  shown  in  Figure  52  contains  several  areas  of  interest  which  may 
create  issues  (tightly  spaced  lines,  sharp  curves).  Primarily  we  expect  there  to  be  a  significant 
increase  in  the  number  of  extracted  segments.  In  addition  to  this  increase  the  completeness 
of  the  extracted  network  may  not  be  as  accurate  as  data  set  one.  This  is  due  to  the  lack  of 
measurements  in  some  areas,  which  causes  the  Hough’s  failure  to  extract  segments  in  the  area. 
The  initial  extraction  is  shown  in  Figure  60. 


Fig.  60:  Data  Set  #2  Phase  One 

The  line  extraction  portion  of  the  algorithm  took  approximately  25  minutes  to  complete.  As 
with  the  first  data  set  the  extracted  line  segment  data  structure  was  stored  before  any  additional 
functions  were  implemented.  Prior  to  any  removal  or  merging  of  line  segments,  the  data  structure 
contains  141  individual  line  segments.  The  next  step  was  the  beginning  of  the  manual  intervention 
to  clean  up  the  road  network.  This  stage  required  the  removal  of  several  line  segments  so  that 
we  could  attempt  to  extract  a  better  fit  for  some  of  the  curved  regions  as  well  as  merging  and 
blending.  Figure  61  shows  the  results  of  the  first  stage  of  manual  interactions  taken. 

The  newly  modified  data  structure  now  contains  105  line  segments  and  57  polynomial  blends. 
With  the  second  phase  completed  the  extraction  was  attempted  again  by  reducing  the  quantization 
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Fig.  61:  Data  Set  #2  Phase  Two 


of  9  from  1  to  0.25  and  the  p  from  2  to  1.  In  addition  the  value  of  the  minimum  length  parameter 
was  reduced  to  10  to  identify  smaller  lines  and  the  fill  gap  threshold  was  reduced  to  30.  The  result 
of  this  extraction  is  shown  in  Figure  62  where  the  discrepancies  can  be  seen  by  comparing  against 
Figure  61,  the  areas  where  the  lines  were  clearly  removed  with  a  large  set  of  measurements  in 
the  left  plot. 
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Fig.  62:  Data  Set  #2  Phase  Three 


An  additional  24  line  segments  have  been  identified  after  the  extraction  process  was 
completed.  Again,  the  user  implemented  a  few  merges  as  well  as  polynomial  blends  in  an 
attempt  to  clean  up  and  more  accurately  define  the  road  network.  This  is  shown  in  Figure  63. 


Fig.  63:  Data  Set  #2  Phase  Four 
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This  newly  modified  data  strueture  eontains  126  line  segments  and  64  polynomial  blends. 
In  the  next  and  final  extraction  phase,  only  the  minimum  length  threshold  was  changed 
from  10  to  5  allowing  for  very  small  line  segments  to  be  identified.  This  resulted  in  the 
identification  of  an  additional  17  line  segments  which  are  much  more  difficult  to  identify  in 
the  plot  due  to  the  small  size.  Figure  64  shows  the  results  of  the  final  extraction  process. 


Fig.  64:  Data  Set  #2  Phase  Five 

Additional  merges  and  blends  were  implemented  and  this  new  structure  contains  139  line 
segments  and  74  polynomials.  Figure  65  shows  the  results  of  the  final  manual  actions  taken. 


Fig.  65:  Data  Set  #2  Phase  Six 


Approved  for  Public  Release;  Distribution  Unlimited. 

103 


A  final  measure  was  taken  due  to  a  few  last  ehanges  in  the  data  structure  and  to  attempt  to 
associate  remaining  data.  The  distance  threshold  was  increased  from  8  pixels  to  10  and  the 
final  results  of  the  manual  interactions  and  extraction  process  are  shown  in  Figure  66. 


Fig.  66:  Data  Set  #2  Final  Extraction  Results 

This  final  data  structure  contains  139  line  segments  and  74  polynomials,  the  same  as  in  the 
previous  phase  since  no  additional  extraction  was  done.  From  stage  one,  the  preliminary  pro¬ 
cessing  of  the  data,  until  stage  seven  approximately  2.5  hours  passed.  Therefore,  to  process  the 
entire  data  set  two  from  start  to  finish  took  approximately  2  hours  and  45  minutes.  A  total  of 
177  measurements  remained  after  all  phases  of  processing  resulting  a  99.8%  of  the  data  being 
associated  with  an  extracted  feature  (88,508/88,685). 

5.3.2  Uncertainty  Analysis 

As  with  the  first  data  set  we  present  here  the  uncertainty  in  the  estimated  parameters  using  the 
various  equations  from  Section  3.11.  First  we  present  the  results  for  the  estimated  CRLB 
matrix  corresponding  to  the  line  estimate’s  parameters.  Figure  67  shows  the  three  sigma  bounded 
region  for  the  lines,  which  again  due  to  the  number  of  measurements  associated  with  each  of 
the  lines,  converges  to  very  finite  values  for  most  lines.  A  few  of  the  lines  have  very  few 
associated  measurements  and  the  bounds  are  very  lax. 

We  supplement  Figure  67  with  an  example  of  a  line’s  CRLB  estimate.  Equation  178  refers  to 
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Fig.  67:  Data  Set  #2  Line  CRLB 
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a  line  whieh  has  1,034  measurements  assoeiated  with  it.  Sinee  there  are  less  measurements  in 
this  data  set  as  well  as  more  lines  when  eompared  against  the  first  data  set,  the  number  of 
measurements  assoeiated  with  a  line  segment  has  signifieantly  deereased  to  around  1,000  at 
most. 


(178) 


We  must  then  use  the  Jaeobian  transformation  to  obtain  the  eovarianee  in  terms  of  x  and  y.  A 
single  measurement’s  eovarianee  matrix  is  given  by  179: 
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The  Jaeobian  transformed  CRLB  is  then  added  to  the  measurement  eovarianee  defined  in  Equa¬ 
tion  179  to  produee  the  final  eovarianee  in  the  measurement  ineluding  the  eovarianee  in  the 
estimate  given  by  Equation  180. 
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Figure  68  shows  a  few  measurement’s  one  sigma  ellipses  with  the  line  estimate. 

One-Sigina  Uncertainty  using  Enor  Ellipses 
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Fig.  68:  Data  Set  #2  Error  Ellipses  One  Sigma 

Using  the  error  ellipses  we  once  again  compute  the  mean  covariance  in  x  and  y.  This  covariance 
is  then  used  to  represent  a  mean  error  in  the  line  using  one  sigma  we  can  see  from  Figure  69 
that  this  encompasses  the  majority  of  the  data  associated  with  each  of  the  line  estimates. 
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Fig.  69:  Data  Set  #2  Mean  Line  Uncertainty 


There  are  no  ellipses  identified  in  this  data  set  but  there  is  a  significant  number  of  polynomial 
blends.  A  polynomial  blend  in  this  data  set  typical  contains  from  400  to  1,000  measurements.  We 
present  the  results  for  a  polynomial  which  contains  1,032  associated  measurements  in  Equation 
181: 
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3.9162e  -  12 
2.9909e  -  10 
-2.0009e  -  08 
-1.5304e-06 


2.9909e  -  10 
2.2842e  -  08 
-1.5281e-06 
-0.0001 


-2.0009e  -  08 
-1.5281e-06 
0.0001 
0.0078 


-1.5304e  -  06 
-0.0001 
0.0078 
0.5980 


(181) 
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6.0  CONCLUSIONS 


Many  concepts  have  been  brought  together  throughout  the  production  of  this  report.  The  Hough 
transform  was  thoroughly  explained  and  derived  in  Section  4.3.4.  In  addition  to  the  standard 
line  Hough  transform  the  circle  and  ellipse  transforms  were  explained  and  the  circle  case  was 
implemented.  The  Least  Squares,  Total  Least  Squares,  and  Least  Squares  ellipse  algorithms  were 
defined  in  Section  4.4.2.  The  purposes  of  each  was  explained  and  the  equations  were  derived 
with  the  appropriate  assumptions.  Section  3.11  derived  the  uncertainty  in  the  estimates 
developed  for  the  line,  third  order  polynomial,  and  the  ellipse.  These  derivations  were 
accompanied  with  examples  and  Monte  Carlo  simulations  which  showed  that  the  solutions 
converged. 

Sections  4.3.4-3.11  developed  concepts  and  equations  which  were  required  during  the 
production  of  the  Graphical  User  Interface  defined  in  Section  4.0.  The  GUI  allows  for 
user  interaction  at  different  stages.  The  accuracy  of  the  extracted  road  network  not  only 
depends  on  the  automatic  extraction  of  the  line  segments  via  the  Hough  transform  but  also  on 
the  time  spent  by  the  user.  As  the  results  showed,  99.98%  of  the  data  in  the  first  data  set  was 
associated  with  geometric  features  and  in  the  second  data  set  the  association  was  99.8%. 

Full  automation  of  an  accurate  road  network  would  be  extremely  complex.  The  ability  to 
identify  an  ellipse  in  an  image  automatically  was  briefly  explored  in  Section  4.3.4,  however  the 
approaches  available  would  require  a  significant  amount  ofcomputation  time  as  they  compare 
every  pixel’s  gradient  direction  or  each  line’s  orientation.  In  addition  to  the  ellipse  extraction, 
the  merging,  blending,  trimming,  and  extending  processes  would  need  to  be  automated.  The 
issues  with  each  of  these  is  that  thresholds  would  need  to  be  established  corresponding  to  the 
line  proximity  and  orientations. 

The  desired  output  of  this  work  was  a  user  interface  which  could  be  easily  utilized  and 
depending  on  the  time  available  would  provide  an  accurate  road  network  (accuracy  directly 
proportional  to  time  spent).  The  results  obtained  as  far  as  associated  measurements  was  very  good 
however  these  results  are  directly  correlated  with  the  amount  of  time  spent  on  the  processing. 
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The  Euclidean  distance  threshold  and  the  endpoint  tolerances  work  well  for  the  two  data  sets 
currently  available.  These  values  may  not  work  well  with  additional  data  sets  and  can  easily  be 
incorporated  as  an  additional  user  input  into  the  GUI  initialization  process. 

The  algorithm  developed  currently  has  no  potential  ability  to  be  computed  in  parallel.  This  is 
due  to  the  recursive  fashion  in  which  a  single  line  is  identified  and  the  points  are  removed 
from  the  image  and  then  another  line  is  identified  with  the  remaining  data.  However  there  are 
areas  in  which  the  algorithm  can  be  cleaned  up  and  better  coded.  In  addition  to  cleaning  up  the 
code,  MatLab  is  not  an  efficient  image  processing  package.  As  seen  with  the  Hough  transform, 
the  MatLab  Hough  function  is  actually  a  compiled  mexfile  (C  code). 

Full  automation  as  stated  previously  would  be  extremely  complex,  however  the  merging 
process  may  be  automated  by  identifying  lines  with  a  similar  orientation  (slope  and  intercept) 
in  addition  to  the  use  of  thresholds. 

Another  possible  area  of  improvement  would  be  with  the  blending  function.  Currently  a  third 
order  polynomial  is  used  which  considers  the  last  xx%,  defined  by  the  user,  of  each  of  the  two 
lines  selected.  Other  options  exists  such  as  higher  or  lower  order  polynomials.  The  third  order 
was  chosen  as  it  appeared  best  for  the  two  data  sets. 

During  the  association  phase  of  the  extraction  process,  the  distance  from  the  current  measure¬ 
ment  to  each  of  the  features  is  computed  and  the  minimal  distance  is  then  found.  This  area  of  the 
algorithm  appears  to  be  sluggish  in  its  implementation  and  further  adaptation  is  required  to 
improve  its  efficiency. 
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List  of  Acronyms 

CHT  Circle  Hough  Transform 
EHT  Ellipse  Hough  Transform 
EaHT  East  Hough  Transform 
EHT  Euzzy  Hough  Transform 
GMTI  Ground  Moving  Target  Indicator 
GPS  Global  Positioning  System 
GUI  Graphieal  User  Interfaee 
ES  Eeast  Squares 
MHT  Multiple  Hypothesis  Tracking 
RMSE  Root  Mean  Square  Error 
SAR  Synthetic  Aperture  Radar 
TES  Total  Eeast  Squares 
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